DNA looping provides stability and robustness to the bacteriophage A switch 
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The bistable gene regulatory switch controlling the transition from lysogeny to lysis in 

bacteriophage A presents a unique challenge to quantitative modeling. Despite extensive 

characterization of this regulatory network, the origin of the extreme stability of the lysogenic 

state remains unclear. We have constructed a stochastic model for this switch. Using 

Forward Flux Sampling simulations, we show that this model predicts an extremely low 

rate of spontaneous prophage induction in a recA mutant, in agreement with experimental 

^ observations. In our model, the DNA loop formed by octamerization of CI boimd to the 
^ ! 

On , Ol and Or operator regions is crucial for stability, allowing the lysogenic state to remain 

. stable even when a large fraction of the total CI is depleted by nonspecific binding to genomic 

CN ■ 

• ■ DNA. DNA looping also ensures that the switch is robust to mutations in the order of the Or 

■ 

' binding sites. Our results suggest that DNA looping can provide a mechanism to maintain 

, a stable lysogenic state in the face of a range of challenges including noisy gene expression, 

nonspecific DNA binding and operator site mutations. 

X 

5t , The bistable developmental switch controlling the transition from lysogeny to lysis in bacterio- 

phage A is one of the best characterized gene regulatory networks In the lysogenic state, the 
phage A genome is integrated into the chromosome of the Escherichia coli host cell and is essen- 
tially dormant, due to expression of the cl repressor gene, the product of which represses cro and 
other genes (Fig. [T]). A transition to the lytic switch state can occur in response to DNA damage 
(via UV irradiation), when CI molecules are cleaved by RecA. Transcription of the cro gene from 
Pr then triggers a cascade of gene activation, leading to phage excision, replication and cell lysis. 
In mutants where this cascade is blocked, a state with elevated Cro levels is stable for several cell 
.e„e,.a.o„. Th. i. .„o„„ a. the a„.-i..„„e .a. Q. A si.p.e i„..i.ive exp>a„a«o„ „as 
been presented for this bistability in the lysogenic state, CI is dominant and cro is repressed; 
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yet, once cro begins to be expressed, Cro dimers repress transcription of c/, making the transition 
to lysis inevitable. However, quantitative measurements have revealed a puzzle: the lysogenic state 
of the phage A switch is both extremely stable and robust to rewiring of its transcriptional 

regulatory interactions Q]. These features have not yet been explained by mathematical models. 
Here, we present dynamical simulations of a stochastic mathematical model that reproduces this 
seemingly mysterious behavior. Our simulations provide evidence that DNA looping plays a key 
role in ensuring the stability and robustness of the switch. 

The molecular interactions controlling the transcription of cl and cro have been studied in 
great detail. CI and Cro bind as dimers to the operator sites O/jl, 0_r2 and 0_r3 (Figured]), which 
control expression of the cl and cro genes from the Prm and Pr promoters. Transcription from 
the unactivated Prm promoter is about 10 times weaker than from Pr; however, when a CI dimer 
is bound at Or2, transcription from Prm is enhanced to about the same level as from Pr [i.e. 
the two promoters can compete with one another only when CI is bound at Or2). CI dimers bind 
preferentially and co-operatively to the OrX and Or2 sites, which overlap the Pr promoter. When 
CI is bound to these two sites, cro {Pr) is repressed and cl [Prm) is activated. Cro dimers bind 
preferentially to 0_r3, which overlaps Prm-, so that when this site is occupied, cl is repressed 
An important additional component of the network architecture is a DNA loop which can form 
between the Or site and the left operator site O^, located 2400bp from Or, which also has three 
adjacent binding sites for CI dimers (OlI) Ol2 and Ol3). The loop is mediated by octamerization 
between pairs of CI dimers bound at Or and Ol Q]- The role of this DNA looping interaction in 
the function of the phage A switch remains a subject of debate [3,1^,0, l3, U]- 

Quantitative measurements have revealed several intriguing features of the phage A switch. 
Firstly, the lysogenic state is extremely stable , despite the stochastic nature of the un- 

derlying gene regulatory network Stochastic fluctuations in gene regulation ("noise") might be 
expected to cause spontaneous transitions from the lysogenic to lytic states, even in lysogens lack- 
ing RecA; yet the rate of these transitions is so low as to be almost unmeasurable. Secondly, recent 
measurements of Prm and Pr promoter activity sugges t that only a small fraction of the total CI 



in the cell is available for binding to Or Q, 1^, 112 
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14( 1 . while other measurements show that the 



total concentration of CI in the lysogen varies dramatically from cell to cell [15|. Taken together, 
these results suggest that the stability of the lysogen is rather insensitive to the number of free 
intracellular CI molecules. Despite this stability, transition to lysis occurs readily in wild-type 
phage on UV irradiation, which leads to cleavage of CI by RecA. Finally, and remarkably, switch 
function is robust to changes in the gene network architecture itself: when the order of the three 
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FIG. 1: Schematic illustration of (a): the On region of the phage A switch, and (b): the DNA loop formed 
by octamerization of CI dimers boimd of Or and Ol- 

Or binding sites is altered so that 0/jl is replaced by 0_r3 or vice versa, the network remains 
functional 

Computer simulations should be an excellent tool for explaining this behavior. Although 
stochastic simulations have successfully been used to model the initial developmental choice be- 
tween lysogeny and lysis for lambda-infected cells modelling spontaneous switching of an 
already established lysogen has proved problematic. Despite the fact that a wealth of biochemical 



parameters are available 
functional robustness 



"or this network, no model has reproduced its extraordinary stability and 
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2l[. In this paper, we present a model that takes into 



account the stochastic character of the chemical reactions and includes DNA looping and depletion 
of free CI and Cro by nonspecific binding to genomic DNA. Our model also explicitly describes the 
detailed dynamics of the binding of transcription factors to the promoters. Accurate computation 
of spontaneous switching rates for this large reaction set is achieved using the Forward Flux Sam- 
pling (FFS) rare event simulation method 



23l |. in combination with temporal coarse-graining 



of dimerization and nonspecific DNA binding reactions. Our simulations show that this stochastic 
model can reproduce the bistability of the switch and its robustness to operator site mutations, 
as well as the extreme stability of the lysogenic state, even in the presence of nonspecific DNA 
binding. 

In this work, we study the effect on switch function of two key parameters: the strengths of 
the DNA looping and nonspecific binding interactions. We find that the DNA looping interaction 
plays a crucial role. In the absence of the looping interaction, a highly stable lysogenic state can 
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be achieved, but this state is very sensitive to depletion of free CI and to operator site mutations. 
When looping is included in the model, the lysogen is insensitive to CI depletion and robust to 
rearrangement of the operator sites. We conclude that DNA looping may play an important role 
in allowing the phage A switch to function reliably even under highly destabilizing conditions in 
the host cell. 



I. THE MODEL 



Our model consists of a set of chemical reactions, simulated using the Gillespie algorithm 
The components of the model are: dimerization of CI and Cro proteins, binding of CI and Cro 
dimers to specific DNA binding sites OrI, Or2, Or3, OlI, Ol2 and Ol3, binding of RNA poly- 
merase (RNAp) to promoters Prm, Pr and Pl, transcription of cl and cro, translation of the 
corresponding mRNA transcripts, degradation of mRNA transcripts and removal of CI and Cro 
monomers and dimers from the cell. Our model also includes formation of a DNA loop between 
Or and Ol, mediated by a CI octamer, and nonspecific binding of CI and Cro dimers to genomic 
DNA. The key parameters that we vary are the strength of the nonspecific DNA binding interac- 
tion AGnsb and the strength of the DNA looping interaction AGioop- Other parameters are fixed 
using biochemical data as far as possible. The model parameters are discussed briefly here and 
described in full in the Supporting Information. 



A. Host cell parameters 



We assume that the E. coli host cell is growing rapidly (doubling time 34 min [3|]), and has 3 



copies of the Or and Ol operators J], in a cell volume of 2/im^ [4|. The concentration of free 



r 



RNAp in the cytoplasm is taken to be 50nM 



25( 1 ■ but our conclusions are not sensitive to this 



parameter, as we demonstrate in the Supporting Information. 



B. Operator binding dynamics 



Equilibrium constants from the literature were used for CI 

mm 

and Cro [28|] binding to OrI, 



Or2, OrZ, for RNAp binding to Prm and Pr 



29( 1 and for CI binding to the Ol sites [8|]. Cro is 



assumed to bind to and Or sites identically. The total number of possible (unlooped) configu- 
rations of the Or and Ol operators are, respectively, 40 and 36. Since we are performing dynamical 
simulations, we require rate constants ka and for association and dissociation. For all association 
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rates, we used the diffusion-limited value ka = AnDa = 0.314/im^s~^ (taking the diffusion constant 
D = 5^m^s~^ and the molecular size a = 5nm). The rate constant for dissociation, in s~^, was 
then deduced from the equilibrium constant, using ka/kd = (exp [— AG/i2T])/(6.023 x 10^)/im^, 
where ka is in ^m^s~^, AG is in kcal/mol, RT = 0.616kcal/mol at 37C, and 6.023 x 10^/Ltm^ is a 
volume conversion factor. 



C. Protein and mRNA production and removal 



We model transcription as a single reaction in which an mRNA molecule is produced when 
RNAp is bound to a promoter. Prm activity is enhanced when a CI dimer is bound at 0_r,2. 
Transcription rates are 0.014s~^ for Pr, O.OOls"^ for unstimulated Prm and O.Olls"^ for stimulated 
Prm 
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3QI, 
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32]. All mRNA transcripts are degraded with a half-life of 2 mins. Translation 
and protein folding are combined into a single step. The model produces a statistical distribution 
for the number of proteins produced per transcript, which is governed by the balance betwen the 
translation and mRNA degradation rates. The average of this distribution (the "burst size") is 6 
and 20 for CI and Cro respectively. The burst size for Cro follows ref. J] . The value for CI is based 
on the observation that the CI ribosome binding site (RBS) is ~ 7-fold weaker than that of LacZ 



33], and that the LacZ burst size is 30-40 
Shean and Gottesman 
where the cell eye 
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351 ]. A somewhat weaker CI RBS was observed by 



half-life 42min 



36]. Protein monomers and dimers are removed with rate constant ln(2)/T, 

n 

e time T =34min [3|]. We also include active degradation of Cro monomers with 



3- 



D. Dimerization 
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38| and — 8.7kcal/mol for Cro 



Dimerization free energies are taken to be — llkcal/mol for CI 
39|. The association reaction is again assumed to be diffusion-limited. To increase the efficiency 
of our simulations, we coarse-grain the monomer-monomer association and dissociation reactions 



for both CI and Cro [40|], as described in the Supporting Information. 



E. DNA looping 

When an Or and an Ol operator each carry at least two adjacent CI dimers (at the 1-2 or 2-3 
sites), these operators can associate to form a "looped state", with rate /cioop, which dissociates 
with rate feunioop- The total number of possible looped states is 49. Binding of CI dimers to the 
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two non-octamerized sites in a loop occurs with a cooperativity factor exp[—AGtct/RT] where 
AGtet = — 3kcal/mol ^j. Because we assume fast growth of the host cell, our model contains 
3 copies of the host genome and 3 copies of the phage A switch. Since the loop is much longer 
than the persistence length of DNA, we assume that any Or — Ol combination can form a loop. 
The strength and dynamics of the DNA loop in vivo are unknown. We therefore test the effects of 
DNA looping on the network behavior, by varying the ratio fcioop/^unioop = exp [—AGioop/ RT]. We 
generally assume a fixed value 62.1s~^ for /cioop (arising from considerations of polymer dynamics 
as discussed in the Supporting Information), but we find that only the ratio is important. 



F. Nonspecific DNA binding dynamics 

We model nonspecific DNA binding [l^ by including in our reaction set association and 
dissociation of CI and Cro dimers to 10^ genomic DNA sites, corresponding to 2-3 copies of the 
bacterial genome. The association rate ka is assumed to be diffusion-limited, and we assume 
identical nonspecific binding affinities for CI and Cro. Nonspecifically bound dimers are removed 
from the cell with rate constant ln(2)/T. We do not model nonspecific binding of RNAp, since 
our value of 50nM corresponds to the free RNAp concentration [25|. To investigate the effects 
of nonspecific DNA binding on the model switch, we vary the parameter AGnsb where ka/kd = 
(exp[-AGNSB/^r])/(6.023 x 10*^)/im3. We assume that these reactions are fast compared to 
the other reactions in the network, and can be coarse-grained [4^ as described in the Supporting 
Information. 



II. BISTABILITY 



Our model represents a mutant phage in which the lytic pathway is nonfunctional, since we do 
not model the lytic genes downstream of cro Such mutants show two stable states: the very 
stable lysogenic state, with high CI and low Cro levels, and the anti-immune state, with elevated 
Cro and little CI, which is less stable, but is nevertheless maintained for several generations Q]. 
For our model to reproduce this bistability, simulations initiated in either of the lysogenic or anti- 
immune states should remain stable, only rarely making a spontaneous transition to the other state. 
Figure [2] shows that this is indeed the case: in Figure [2^, a simulation initiated in the lysogenic 
state remains in that state, while in Figure [2}3, a simulation run with the same parameters but 
initiated with a high Cro concentration and little CI remains in the anti-immune state. Figure [3] 
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FIG. 2: The model switch is bistable. Typical simulation trajectories are shown for AGioop = — 3.7kcal/mol 
and AGnsb = — 4.1kcal/mol. (a): A simulation initiated with 150nM CI and no Cro molecules remains in 
the CTrich state for many hours (b): A simulation initiated with 400nM Cro and no CI molecules remains 
in the Cro-rich state. 



shows the range of parameter values for which our model shows bistability. Our simulations give 
steady state concentrations of CI and Cro in good agreement with measured values for the lysogenic 



and anti-immune states, 
a measured value of 220 



or the lysogenic state, we obtain ~ 200 — 400 CI per cell, compared to 



4l[ |. For the anti-immune state, Cro per cell ranges from 250 — 900 



(depending on AGnsb); corresponding to concentrations of 200 — 750nM, compared to a measured 



concentration of ~ 400nM [1^. These values are presented as functions of AGioop and AGnsb in 
the Supporting Information. In our model, the DNA looping interaction decreases the lysogenic 
CI concentration by as much as a factor of 2, in agreement with the observations of Dodd et al 
0, We have also simulated mutants without the Or3 or the 0x3 binding sites, corresponding 
approximately to the 0R3-rl and OL3-4 mutants of refs Q], Q] 



and [ll[. For these mutants, we 
find lysogenic CI concentrations 1.9 and 1.8 times the wild-type values (for AGioon = — 3.7kcal/mol 
and AGnsb = — 2.8kcal/mol), in reasonable agreement with the results of ref [8| (factors of 2.8 



and 3 respective^ 
identified in ref 



y), even though our model does not include the up-regulation of Prm by the loop 



111 ], or the effect on Prm of the 0R3-rl mutation. For the same parameters, the 



Pr promoter is repressed by a factor of 1.6 in the anti-immune state, in good agreement with ref. 



III. EXTREME STABILITY OF THE LYSOGENIC STATE 



The spontaneous switching rate from lysogeny to lysis in recA- mutants is so low as to be 



almost undetectable, being less than 2 x 10 ^ per cell per generation 



4]- Reproducing such 



low spontaneous switching rates while maintaining bistability provides an extreme challenge for a 
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FIG. 3: Bistability of the model switch as a function of the nonspecific binding strength AGnsb and the 
DNA looping strength AGioop- Blue squares represent parameter sets for which only the lysogenic state 
is stable (no stable anti-immune state); for red diamonds only the anti-immmie state is stable (no stable 
lysogen), while green circles represent parameter combinations where the model switch is bistable. 



a 



computational model [J]. We quantified the stability of the lysogenic and anti- immune states for 
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231]. This 



our model, using the forward flux sampling (FFS) rare event simulation method 
method allows us to compute the rate of spontaneous switch flips, even though such flips would 
hardly ever be observed in a typical "brute-force" simulation run. FFS uses a series of interfaces 
(defined by an order parameter) between the initial (Cl-rich) and final (Cro-rich) states to split 
a switching event into a number of more likely transitions between successive interfaces. Here, 



our order parameter was the difference between the total number of Cro and CI molecules^ 
method has previously been used for a simple model of a mutually repressing genetic switch 
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'his 



231]. 



Details of the FFS method and its implementation are given in the Supporting Information. 

Our results are listed in Table [H In the absence of DNA looping (top two rows of Table H]), 
the lysogen is stable for ~ 10^ generations only when nonspecific DNA binding is absent. When 
nonspecific binding is included in the model, the lysogen becomes much less stable: a relatively 
modest nonspecific binding strength AGnsb = — 2.1kcal/mol produces a lysogen that is only stable 
for 1700 generations: approximately a million times less stable than the experimental lower 
bound. In contrast, when DNA looping is included in the model (bottom three rows of Table 
H]), extremely stable lysogens can be achieved for a wide range of DNA looping and nonspecific 
DNA binding strengths. These lysogenic states are in fact even more stable than those observed 
experimentally Q]. One possible explanation might be the effects of passing DNA replication forks 
(not included in the model), which might be expected to destabilize looped DNA configurations 
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and/or remove bound CI from the operator sites. The anti-immune state for our model is much 
less stable than the lysogenic state, in agreement with experimental observations P, . 

TABLE I: Spontaneous switching times (inverse of calculated switching rates) from the lysogenic to anti- 
immune states and from the anti-immune to lysogenic states, for wild-type phage A, computed using FFS. 



AGloop 


AGnsb 


Switching time 


Switching time 


kcal/mol 


kcal / mol 


lysogen —>■ anti-immune 


anti-immune — > lysogen 






(generations) 


(generations) 


no loop 


no NSB 


(3.6 ±0.1) X 10^ 


2300 ± 100 


no loop 


-2.8 


1700 ± 50 


(7± 1) X 10^ 


-5.2 


-4.1 


< 1023 


43 ± 1 


-3.7 


-2.8 


(3 ± 2) X 10^5 


26 ± 1 


-1.0 


no NSB 


(6 ± 2) X 10^'^ 


130 ± 10 " 



"A cell generation time of 34 min is assumed. FFS calculations used 8-85 interfaces, 1000-10000 configurations at 
the first interface and 500-10000 trials per interface, and were averaged over 10-100 runs. 



IV. DNA LOOPING ALLOWS THE SWITCH TO FUNCTION DESPITE CI 

DEPLETION 

It is believed that a large fraction of the transcription factors in a bacterial cell are unavailable 
'or binding to their specific binding sites because they are nonspecifically bound to genomic DNA 



42l |. Recent expression measurements for the Prm promoter have suggested that this is the 
case for CI This nonspecific binding poses a severe challenge to the phage A switch. 

Although both CI and Cro are depleted by nonspecific binding, the Prm promoter is intrinsically 
weak, and requires activation by CI at Or2 to compete effectively with Pp. One would therefore 
expect nonspecific DNA binding by CI to drastically destabilize the lysogen, and to compromise the 
bistability of the switch. Table HIl of the Supporting Information shows how the concentration of free 
CI depends on the nonspecific DNA binding strength. To characterise the effects on switch function, 
we investigated the range of nonspecific binding strengths over which our model gave bistability, 
for different values of the DNA looping parameter AGioop- Our results are shown in Figure [3l 
In the absence of DNA looping (top row of Figure [3]) , bistability is indeed strongly compromised 
as the parameter AGnsb is increased in magnitude (left to right). For | AGnsb I >~ 3kcal/mol, 
the model is no longer bistable; the lysogenic state cannot be sustained and only the anti-immune 
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state is stable. However, when the DNA looping interaction is included in the model, bistability 
is maintained over a much wider range of AGnsb values. In figure [3l the width of the green 
(bistable) region increases dramatically as the parameter |AGioop| increases. Our model therefore 
suggests that one role of the DNA looping interaction may be to ensure that the switch continues 
to function even when the level of intracellular free CI is depleted by nonspecific DNA binding 
0, ^ or by cell-to-cell fluctuations 15]. We note that this stability to CI fluctuations is not 
expected to prevent lysis from occurring on UV irradiation of the wild-type phage: even for a 
strong looping interaction, rapid degradation of CI by RecA will eventually lead to so little CI 
being present that the loop cannot be maintained, upon which lysis will occur. To check this, 
we simulated a version of the model in which we fixed the total CI concentration, for a typical 
parameter set (AGioop = — 3.7kcal/mol, AGnsb = — 2.8kcal/mol). When we artificially lower the 
CI concentration to ~ 10% of the lysogenic steady-state level (30 CI molecules per cell), the system 
flips to the anti-immune state. This result is in good agreement with the observation of Bailone et 
al 43] that prophage induction occurs at a CI concentration about 10% of that in the lysogen. 



V. DNA LOOPING CAUSES ROBUSTNESS TO OPERATOR MUTATIONS 

In an important series of experiments, Little et al showed that the basic functions of the phage A 
regulatory network are robust to changes in network architecture. Mutants 0/j(121) and Oij(323), 
in which the order of the O^jl, Oi?2 and 0_r3 binding sites was altered compared to the wild- type 
Oij(321), formed stable lysogens (although less stable than the wild-type) which could be induced 



to enter the lytic pathway on UV irradiation 

El 



a. 



To our knowledge, this robustness has not been 

reproduced in computer models [4|. 

We tested whether our model was able to produce stable lysogens for the 0^(121) and 0/j(323) 
mutants, which were created in our simulations by changing the operator binding site affinities. 
We neglect possible changes in the properties of Prm whose DNA sequence is also affected by 
the Little et al substitutions. The range of parameters for which our model mutants are bistable 
is shown in Figure SI We find that the DNA looping interaction plays a key role in ensuring 
robustness. In the absence of DNA looping (top row of Figure H]), our model produces a stable 
lysogen for 0^(121) only if the nonspecific binding strength AGnsb is below a critical value, and 
cannot sustain a stable lysogen for 0^(323) at all. However, when DNA looping is included in the 
model (lower rows of Figure H]), the Or{2>22>) mutant can achieve a stable lysogenic state, and the 
Or{\2\) lysogen becomes able to tolerate stronger nonspecific DNA binding interactions. In fact. 



11 



Op(i2i; 



Op(323) 



no loop- • • • 
0- 

: ■ ■ • • • 

-4^ ■ 



■ • •• 



noNSB 2 3 4 



1 

■ ♦ 


1 1 1 1 1 
♦ ♦ ♦♦ 


: • 


♦ ♦♦♦ : 


: • 


• ♦ ♦♦ : 


L • 


• • ♦♦ J 


1 


• • 

" 



o 
<1 



noNSB 2 3 4 



FIG. 4: Range of bistability of the model switch as a function of the nonspecific DNA binding parameters 
AGnsb and AGioop, in the presence and absence of the DNA looping interaction, for the Little mutants. 
Symbols have the same meaning as in Figure [3l 

if the looping strength is too strong, our model predicts that the Or(121) mutant may lose the 
ability to sustain a stable anti-immune state. 

The steady-state concentrations of CI in the lysogen predicted for these mutants are in rea- 
sonable agreement with the observations of Little et al. For a typical parameter set (AGioop = 
— 3.7kcal/mol, AGnsb = — 2.8kcal/mol), we obtained steady-state lysogenic CI concentrations 
38% and 81% of the wild-type values for Or{\2\) and Oi?(323) respectively, in comparison to 
25-30% and 60-75% measured by Little et al. Calculation of the spontaneous switching rates for 
the 0/j(121) and Or(323) mutants, using FFS, also shows results qualitatively in agreement with 
those of Little ai [3| (see table of results in the Supporting Information). For all combinations of 
AGioop and AGnsb, the 0/j(323) lysogen is less stable than Oij,(121) which is in turn less stable 
than the wild-type, although quantitatively the magnitude of this destabilization is greater than 
the factors of 10 observed by Little et al. Our model allows us to make the tentative prediction that 
(in a recA- host) the stability of the anti-immune states will be 0/j(323) > wild-type > Or{121). 
These stabilities have not yet been measured 



VI. ALTERNATIVES TO LOOPING 



Our results indicate that the DNA looping interaction allows the phage A switch to maintain 
an extremely stable lysogenic state, which is robust to operator mutations, while retaining its 
essential bistability. Could this result have been achieved by evolution in any other way? Recent 
work by Babic and Little 481] has shown that a triple mutant lacking CI cooperativity but with 



increased Prm activity and enhanced binding of CI to 0^2, can maintain a stable lysogen, with 
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increased CI levels compared to the wild-type, and can switch into the lytic state. It is very 
unlikely that this mutant can form a DNA loop [45|. This suggests that increased lysogenic CI 



levels could substitute for cooperative interactions (including DNA looping). To test whether our 
model could reproduce these experiments, we modified our simulation parameters to represent the 
two triple mutants A cl Y210N prm252 Or2up and A cl Y210N prmup-1 OrIup. In both cases, 
we removed all cooperative interactions between CI dimers, including looping, and increased the 
affinity of CI for 0/j2 by a factor of 5. For the prm.252 mutation, the basal and stimulated Prm 
transcription rates were increased by factors of 6 and 2.5 respectively, and for the prmup-1 mutation 
these factors were 10 and 5 48(. With nonspecific binding strength AGnsb = — 2.1kcal/mol, both 
these "mutants" formed stable lysogens, with CI levels ~ 3 and 6 times that of the "wild- 
type". The mutant representing A cl Y210N prmup-1 Or2up also formed a stable lysogen for 
^Gnsb = — 2.8kcal/mol, but A cl Y210N prm252 Or2up did not. These results are in reasonable 



48l |. Furthermore, our results support the view 



agreement with those reported by Babic and Little 
that increased lysogenic CI levels could provide an alternative mechanism for maintaining a stable 
lysogenic state, although this mechanism has apparently not been selected by evolution. 



VII. DISCUSSION 



The phage A switch represents an important test for computational modeling of gene regulatory 
networks. This network has been a paradigm in molecular biology for many years: its architecture 
and biochemical parameters have been extensively studied and its behavior has been thoroughly and 
quantitatively characterized. Nevertheless, computational models have failed to produce results in 
agreement with experimental observations. If modeling is to prove a useful tool in the analysis of 
larger and more complex gene regulatory networks, it is essential that it should produce convincing 
results for phage A. 

The computer simulation model presented here reproduces for the first time both the bistability 
of the phage A switch and its extremely low spontaneous flipping rate. Our model includes only 
known biochemical interactions and uses as far as possible biochemically measured parameters. 
The forward flux sampling (FFS) rare event sampling method allows us to quantify the stability 
of the lysogenic and anti-immune states, even though spontaneous switching rates would never be 
observed in a typical brute-force simulation run. An important difference between our model and 
previous theoretical studies is that we simulate the full transcription factor-DNA binding dynamics. 
These "operator state fluctuations" were eliminated in previous models, using the physical-chemical 
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29( 1 ■ Our previous work 



49|, as well as 



quasi-equilibrium assumption of Shea and Ackers 

that of Walczak and Wolynes 3], has shown that operator state fluctuations can drastically 
affect the switching pathways and the switching rate for bistable genetic networks. Moreover, 



Vilar and Leibler 



47l | have shown that operator state fluctuations are coupled to DNA looping 



dynamics. We therefore model explicitly both the protein-DNA binding dynamics and the DNA 
looping dynamics. Our model also includes nonspecific DNA binding, as does ref. J]. Although 
DNA looping and nonspecific binding are modeled in a simplified manner, our model involves 
many reactions and is computationally expensive to simulate. This problem is alleviated by the 
FFS method, in combination with the coarse-graining of dimerization and non-specific binding 
reactions. We note that a semi-analytical approach, based on Large Deviation Theory, which has 
been applied successfully to a simplified model of the phage A switch {lo| , may provide a promising 
alternative to direct simulations as performed here. 

Our results suggest a key role for the DNA looping interaction in ensuring that the network 
retains its bistability even when exposed to perturbations. Our model requires DNA looping to 
achieve bistability in the presence of nonspecific DNA binding and/or operator state mutations. 
The Prm promoter is intrinsically weak and requires CI to be bound at Or2 in order to compete 
effectively with Pr. The DNA loop enhances CI occupancy at Or2, providing a way to achieve 
sustained activation of the Prri promoter, even when the free intracellular CI concentration is 
very small. We also find, in agreement with Dodd et a/,jthat the presence of the loop increases 
autorepression of Prm by enhancing CI binding at Or3 12]. Looping thus increases the strength 
of both auto-activation via Or2 and auto-repression via OrS; this reduces the fiuctuations in CI, 



50(]. Our model does not include the loop-mediated 



which enhances the stability of the switch 
increase in maximal Prm activity recently discovered by Anderson and Yang 
including this in the model would only strengthen our conclusions. 



ll|; we expect that 



48( 1 show that a stable lysogen can be 



Our simulations and the experiments of Babic and Little 
obtained in the absence of DNA looping and other cooperative interactions, by raising the level of 
CI. This observation is perhaps not so surprising, since the stability of genetic switches is predicted 



5C 



52] . What is then 



to depend exponentially on the expression levels of gene regulatory proteins 
the role of DNA looping? DNA looping not only reduces fluctuations in CI levels, as discussed 
above, but also reduces operator state fiuctuations, which have been predicted to limit the stability 
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of genetic switches 
to be achieved with lower CI leve 
to the host cell of producing CI 



49( 1 . This suggests that the looping interaction allows a stable lysogen 
s. One advantage of this may be a reduction in the energetic cost 



51 



53l ]. Alternatively, the dynamical pathways for the transition 
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to lysis may be more favorable for a switch with DNA looping. 

Several predictions emerge from our simulations. Firstly, we predict that a mutant phage which 
cannot form the DNA loop will form a lysogen with much reduced stability compared to the wild- 
type, and may not be able to sustain a lysogen at all. We further predict that depleting free CI 
(for example by introducing a large number of specific binding sites) will strongly destabilize the 
lysogen in a non-looping mutant but will have a much less severe effect on the wild-type lysogen. 
Finally, our simulations allow us to tentatively suggest an order for the stability of the anti-immune 
states in a non-lysing version (for example Ofj,{121)N~0~) of the Little mutants in a recA- host. 

This work presents an encouraging picture of the potential for computational modeling to 
unravel the contribution of different biochemical mechanisms to observed biological behavior. Using 
a simplified representation of the core part of the phage A switch, we are able to obtain behavior 
in qualitative and partial quantitative agreement with experimental results. Future work should 
include more sophisticated models for DNA looping and nonspecific binding, explicit modeling of 
cell growth, DNA replication and cell division, as well as detailed characterization of the switching 
pathways. Such work should make stochastic modeling, in combination with experiments, an 
important tool in unraveling the mechanisms behind the complex behavior of biochemical networks. 
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VIII. SUPPLEMENTARY MATERIAL 

In section IIXI we discuss in more detail the technical aspects of our stochastic simulations. In 
section IXl we present results including steady-state protein concentrations, depletion of free CI by 
nonspecific binding, promoter activity as a function of CI concentration, switching rate calculations 
for the Little mutants and effects of increasing the RNAp concentration. In section IXIl we test the 
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sensitivity of our conclusions to variations in the parameters of the model. Finally, in Appendix 
IXIH we give a complete list of the chemical components and reactions in our model. We note 
that many of these reactions are permutations of each other, so that the number of rate constant 
parameters is much smaller than the number of reactions. 

IX. SIMULATION DETAILS 

In section IIX Al we discuss how DNA looping is represented in our model. In sections IIXBI 
and IIX CI we give details of how we coarse-grain the dimerization reactions and the nonspecific 
DNA binding reactions in the model. In section IIXDI we discuss Forward Flux Sampling and its 
application to this model. 

A. Modeling DNA looping 

Representation of looping in our simulation model 

In our model, the Ol operator has three binding sites, O^l, Ol'^ and 0^3, which can bind 
CI or Cro dimers, as well as a promoter, Pl, which overlaps the binding site O^l, and which can 
bind RNA polymerase. It is important to include RNAp binding to Pl, because this excludes CI 
binding to 0^1 (even though no gene is expressed from Pl in the model). Binding constants for 
Cro dimers to the three Ol binding sites are assumed to be the same as for Or. For CI, we use 
the binding constants measured by Dodd et al [s]. 

In our model, there are 3 copies of each of the Or and Ol operators in the cell (since we assume 
a doubling time of 34 minutes) . We assume that any of the copies of O r can form a loop with any 
copy of Ol, since the 2400bp separation between Or and Ol on the A genome is long compared 
with the persistence length of DNA (~ 50nm = 150bp). The system can form a loop if both the 
Or and Ol operators are bound by at least two adjacent pairs of CI dimers {i.e. in the 1-2 or 
2-3 positions). Fully bound operators (with dimers in all 3 positions) can also form loops. We 
classify looped states into four categories, according to which binding sites are involved in octamer 
formation. Loops where the octamer is formed by CI dimers at OrI, Or2, Ol^ and Ol2 are 
denoted OLRl, loops whose octamer is formed by dimers at OrI, Or2, Ol2 and Ol3 are denoted 
0LR2, loops whose octamer is formed by dimers at Or2, Or3, OlI and Ol2 are denoted OLR3 
finally loops with octamers formed by CI dimers at Or2, Or3, Ol2 and 0^3 are denoted OLRA. 
For each of these categories, the additional two binding sites, not involved in octamer formation. 
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may also be occupied by CI or Cro dimers (or in some cases RNAp). We therefore designate any 
particular looped state by OLRX{YZ), where X is 1, 2, 3 or 4, as above, and Y and Z denote the 
species that are bound to the two non-octamerized sites on Ol and Oji respectively. For example, 
a looped state where CI dimers bound at 0_r1, Or2, OlI and 0^2 participate in an octamer, and 
additional CI dimers are bound at 0^3 and 0^3 would be denoted as OLRl{RR) in the reaction 
scheme in Appendix IXIII of the Supporting Information (here, "R" denotes "CI repressor"). 

Loop formation is represented in our reaction scheme (see Appendix IXIip by reactions in which 
unlooped states with Cl-bound Or and Ol operators associate to and dissociate from the appro- 
priate OLR states. Association to a looped state always occurs with rate /sioopi and dissociation 
always occurs with rate /cunioopj regardless of which combination of CI dimers and/or RNAp is 
bound. If a looped state has vacant binding sites, binding to these is assumed to be diffusion- 
limited. To take account of the cooperative binding interaction due to tetramer formation between 
the two non-octamerized CI dimers in the fully occupied looped state Q], we reduce the dissociation 
rate of a single dimer from the fully occupied looped states 0LR1{RR) and OLRA[RR) by a factor 

. This scheme for including the cooperativity due 



exp [—lS.Gici/ RT] where AGtet = — 3kcal/mol 
to tetramerization ignores the possibility that loops may form in either orientation 
believe that this simplification is likely to have only a minor effect on our results. 



Ill ], but we 



Estimation of DNA looping parameters from polymer dynamics 

The strength of the DNA loop between Or and Ol in vivo is unknown. We therefore vary the 
ratio /cioop/^unioop = exp [—AGioop/RT] in our simulations, maintaining a fixed value of 62.1s~^ for 
/cioop and varying feunioop [as for protein-DNA binding, we assume that the binding affinity affects 
only the "off" -rate]. In some cases, where this would result in extremely slow loop dissociation 
(slower than the cell doubling time), we increase /cioop to prevent the looping dynamics becoming 
very slow. In all the cases we have examined, we find that only the ratio A;ioop/^unioop is important, 
indicating that in our model, DNA looping/unlooping is not the rate limiting step for switch 
flipping. 

Although we use AGioop, and hence /cunioopi as an adjustable parameter in our simulations, we 
can gain some insight into likely values for /cioop from a consideration of polymer dynamics. For 
simplicity, we assume that the association time for two Cl-bound operators, separated by 2400bp, 
to form a loop, is of the same order of magnitude as the mean relaxation time of a 2400bp DNA 
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molecule, whose ends have been brought into contact. This is given by Meiners as 



(1) 



TTkhTln{LQ/d) 

In this expression, r] is the viscosity coefficient for the cytoplasm, Lq = 24006p = SOOnm is the 
contour length of the DNA, Ip is the persistence length and d = 2.5nm is the thickness of DNA. 



Assuming the viscosity of the cytoplasm to be ten times that of water [55l ]. ~ 10 x 10"'^ Pa.s, 
we obtain ta = 0.016s. This leads to ^loop = ~ 62.1s, which is the typical value used in our 
simulations. We then vary fcunioop according to the chosen value of AGioop- 

In our simulations, we vary AGioop over a range larger than the value of — 0.5kcal/mol obtained 
by Dodd et al by fitting promoter activity data for Prm and Pr. We assess the extent to which 
our model fits this promoter activity data in Section [X CI As discussed by Dodd et a/ [8!], the true 
value of AGioop in vivo is unknown. In vitro measurements of the free energy of association of two 
CI tetramers in the absence of DNA give a value of — 9.1kcal/mol [5(|, but this may be different in 
the presence of DNA, and in any case must be offset by a contribution due to the entropy of DNA 
looping. 



B. Coarse-graining dimerization 

The monomer-dimer association and dissociation reactions: 

CI + CI ^Ch (2) 
Cro + Cro Cro2 

are responsible for the vast majority (~ 99%) of the total computational effort when our model 
switch is simulated in the absence of looping and nonspecific DNA binding. This is because 
these reactions have very high propensities, due to the large number of Cl/Cro molecules in the 
lysogenic/anti-immune stable states. We can greatly increase the efficiency of our simulations, 
without jeopardizing their accuracy, by "coarse-graining" these reactions. We assume that the 
monomer-dimer association and dissociation reactions are fast enough that they reach steady state 
on the timescale of the other, slower reactions in the system (although there is some evidence that 
Cro dimerization kinetics may be slow js^). We can then remove these reactions from our reaction 



set 



as in our previous work on a simpler genetic switch model 



40l ? ]. To do this, we 



18 



define new simulation variables: 

"ci = '^ci + 2nci2 (3) 

"'Cro = "Cro + 2ncro2 

where nci represents the number of CI monomers, ncij the number of CI dimers and the sum of 
these numbers (and likewise for Cro). We rewrite our reaction scheme in terms of the new "chemical 
species" CI and Cro, whose numbers remain unchanged by the association/dissociation reactions 
([2]). Translation reactions now produce CI or Cro, rather than CI or Cro, but with propensity 
computed as before. Protein removal from the system by dilution now removes CI or Cro, also with 
reaction propensity unchanged. The active degradation of Cro monomers now becomes degradation 
of Cro. Here, the new reaction propensity is given by the old rate constant, multiplied by the 
average number (racro)nc-o of Cro monomers that would be obtained by a simulation of the fast 
association/dissociation reactions, at fixed n^-^ [in what follows, we will use angular brackets with 
a subscript to denote an average value of the quantity in the brackets, taken over a simulation 
in which the subscripted quantity is held constant]. All dimer-DNA dissociation reactions now 
produce CI or Cro instead of CI2 or Cro2, with propensity computed as in the full reaction scheme. 
All dimer-DNA association reactions are now represented by the association of two units of CI or 
Cro to Oij or O/,, with propensities equal to the original reaction rate, multiplied by the average 
number {nci2)n^^ or (?^Cro2)r^.|-,-^J of dimers given by the association/dissociation reaction set ([2]) 
for fixed n^jj or Hq-^. To implement this scheme, we need to know (nci)^,^^, ("Cro)^^^-^! ('^Cb)^^; 
and {ncro2)n^-^^, for all possible values of CI and Cro. These values are obtained prior to the main 
simulation run, by numerical solution of the chemical master equations corresponding to Eq.([2]), 
for fixed CI and Cro 40j]. The results are stored in a table for all values of CI and Cro. This table is 
then used to compute the necessary reaction propensities in our simulations of the coarse-grained 
phage A reaction set. 

We have demonstrated this procedure [5^ . Is^ in previous work on a model bistable switch 
formed of two mutually repressing genes [40|. In that work, we showed that coarse-graining the 
monomer-dimer association/dissociation reactions had little effect on the switch flipping rate - in 
contrast to the dimer-DNA association/dissociation reactions, which could not safely be coarse- 
grained. We have also verified that for our phage A model (in the absence of DNA looping or 
nonspecific binding, for simplicity), the lysogenic to anti-immune switching rate computed with 
explicit monomer-dimer association/dissociation is indistinguishable from the same rate as com- 
puted with the coarse-grained reaction set. 
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C. Coarse- graining nonspecific DNA binding 

Nonspecific binding of CI and Cro dimers to E. coli genomic DNA is included in our model via 
a set of association-dissociation reactions to genomic DNA sites ( "D" ) , of which we assume there 
are 10'' (about 70% of the total base pairs in the three copies of the E. coli genome that we suppose 
to be present). These reactions are: 

CI2 + D ^ DCI2 (4) 
Cro2 + D ^ DCro2 

These reactions do not necessarily represent the actual cellular process which results in depletion 
of free intracellular CI and Cro. Such processes are likely to be very complex: these reactions 
simply provide a convenient and simple way to deplete the available CI in our simulations. The 
nonspecifically bound species DCI2 and DCro2 are assumed to be subject to dilution due to cell 
growth, which we model by removing them from the system with rate constant ln(2)/T, where 
T = 34 min is the cell cycle time. The free energy AGnsb of binding of a CI or Cro dimers to a 
D site is varied systematically in our simulations. We assume this parameter takes the same value 
for CI and Cro. The association rate ka is assumed to be diffusion-limited, and the dissociation 
rate is calculated as described in the main text. However, due to the large number of D sites, 
the propensities for these reactions are very much larger than those for any other reactions in the 
system. Direct Gillespie simulation of these reactions is not computationally feasible. Instead, we 
use a coarse-graining scheme which we construct according to the principles described above for the 
dimerization reactions, and in reference [4^. In this scheme, we coarse-grain both the dimerization 
and nonspecific binding reactions. This implies that we assume that both association/dissociation 
to D sites, and association/dissociation of dimers, occur on much faster timescales than the other 
reactions in the system. 

To carry out this coarse-graining, we define new simulation variables: 



'^cr = nci + 2nci2 + 2nDci2 (5) 

"■Cro' ^ "■Cro + 2nCro2 + 2nDCro2 

which represent the total number of CI and Cro molecules in the system, excluding only those 
which are specifically bound to the operators Or and Ol. Our reaction scheme in fact remains 
exactly the same as that described in section IIXBI above, except that CI and Cro are replaced 
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by CI' and Cro', and the averages {nci)n^^, ("-Cro)^^-^, {nch)n^^ and (ncroa)?!^-^ are replaced by 
("-ci)nciM (f^Cro)ncro" ("Ci2)nci/ and (ncro^) n^,, ■ These latter averages must be computed over 
both the dimerization reaction set ([2]) and the nonspecific binding reaction set (j4]), for fixed values 
of Hqii and ?^cro'• We assume that the reaction sets for CI and Cro do not couple {i.e. that the 
number of D sites is very large), so that we can compute the two sets of averages independently. 
We carry out a preliminary set of simulations, in which only reactions ([2]) and (jH) are simulated, 
for either CI or Cro, for fixed or riQ^.^/. Using these simulations, we evaluate the necessary 
averages, which we tabulate for all values of n^j/ and ?^cro'■ We then use these tabulated values 
to compute the required propensities in our coarse-grained simulations of the phage A switch. 
This coarse-graining leads to a speedup of our simulations by a factor of at least 100. Without 
this coarse-graining, the computations described here would have been beyond the limits of our 
computational resources. 



D. Forward flux sampling for the phage A switch 

Forward Flux Sampling (FFS) allows the calculation of rate constants and transition paths for 
rare events in equilibrium and nonequilibrium stochastic simulations. FFS is described in detail 
in refs 



231 and 



23( 1 . It will be sketched briefly here in the context of the transition from the 
lysogenic to anti-immune states of our model. We define a collective variable of the system, or 
order parameter, that distinguishes the initial and final states. We do not assume that this order 
parameter is the true "reaction coordinate" for the transition. For the lysogenic to anti-immune 
transition, our order parameter N is taken to be the difference between the total numbers of CI 
and Cro molecules in the system: 

= total number of Cro — total number of CI (6) 

In the anti-immune state, Cro dominates and N is positive. In the lysogen, CI dominates and N 
is negative. We define values A'^iys and A^'ai such that if < A^iys, the system is assumed to be in 
the lysogenic stable state, and if > A'^aii the system is assumed to be in the anti-immune stable 
state. The chosen values of Niys and A'^a.i vary for different parameter values, since the numbers of 
CI and Cro in the two stable states are affected by the choice of parameters (as shown in Figure 
[5]). The value of Niys should be such that, when the system is in the steady state, fluctuations 
quite frequently take it into the region where N > Niy^. 

We note that although we have chosen N as our order parameter, we could have made an 



21 



Lysogenic 



Anti-immune 




J I I L 

-2 no loop 




200 ■ 



-2 



B 

L 

no loop 



FIG. 5: Steady-state total number of molecules of CI in the lysogenic state and of Cro in the anti-immune 
state, as a function of the DNA looping strength AGioop, for different values of the nonspecific binding 
strength AGnsb- The lysogenic CI concentration is in agreement with measured values and is rather 
insensitive to nonspecific binding, although DNA looping causes a decrease in the lysogenic CI concentration. 
The concentration of Cro in the anti-immune state is insensitive to DNA looping strength (as expected) but 
increases with the strength of the nonspecific binding interaction. 



alternative choice - for example, the total number of CI molecules, or some weighted combination 
of the numbers of CI and Cro. Providing these other choices give good separation of the initial 
and final states, the same rate constants and transition paths should be obtained. However, 
the computational efficiency and sampling effectiveness may be affected by the choice of order 
parameter. 

The aim of the computation is to compute the rate k at which the system makes transitions 
from the region where < A'^iys to the region where N > A^ai- This is the frequency with which 
simulation trajectories leave the lysogenic state and enter the anti-immune state. This frequency 



is very low, but it can be re-expressed as 



k = $ivsPi 



lys^lys— +ai 



(7) 



Here, <l>iys is the flux of trajectories out of the lysogenic state, or the frequency with which the 
system crosses the interface N = Niys coming from the lysogenic state. This flux must be mul- 
tiplied by the probability Piys_>ai that one of these trajectories which leaves the lysogenic state 
will subsequently reach the anti-immune state rather than returning to the lysogenic state. The 
flux $iys can be calculated very accurately because the system makes many crossings of = A^iys- 
The probability Piys^ai is however very small and difficult to calculate. This problem can be over- 
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come by placing a number of "interfaces" between the lysogenic and lytic states, at values of N 
intermediate between A^iys and A^'ai- We define A^'o = iViyg, = N^\, and place interfaces Ni for 



1 < i < n — 1 between A^'o and Nn- Eq. ([7]) can then be rewritten as 

n-l 

k = <^iy,l[P{Ni+i\Ni) (8) 

i=0 

where P{Ni^i\Ni) is the probability that a trajectory which has reached interface Ni, coming from 
the lysogenic state, will subsequently reach the next interface Ni^i, rather than returning to the 
lysogenic state. Since the interfaces may be placed arbitrarily close together, the probabilities 
P{Ni-^i\Ni) can be large, and thus easy to compute accurately. Once these probabilities and the 
flux $iys have been computed, the rate for the lysogenic to anti-immune transition can be obtained 
using Eq.dSI). 

We compute the flux <I>iys using a simulation run initiated in the lysogenic state. Every time 
the simulation crosses = A'^iys in the direction of increasing N, we increment a counter. At the 
same time, we store the configuration of the system at the moment that the crossing occurs. At 
the end of this simulation run, we divide the value of our counter by the total simulation time to 
obtain 'I'lys- We have also obtained a collection of system configurations corresponding to crossings 
of AT = A^iy,. 

We now turn to the computation of the probabilities P{Ni^i\Ni). To compute P(A'i|A'^o)i we 
fire many "trial runs" from the collection of configurations that we have obtained at Nq. In each 
trial run, we choose a configuration at random from our collection and use it as the starting point 
for a new simulation run, which is continued until either the system reaches A'^i, or it returns to 
the lysogenic state A^ < Nq. If A^i is reached, the final configuration of the trial run is stored in a 
new collection. We repeat this trial run procedure many times, and finally estimate P(A^i|A'^o) as 
the fraction of trial runs that successfully reached A^i . We now use the collection of configurations 
which we have just collected at A^i to compute P(A'^2|A^i): once again we fire many trial runs 
which are continued until either A'2 or A'o is reached. We repeat this procedure for each interface, 
until we have computed all the P{Ni^i\Ni) values. At this point, trajectories corresponding to 
transitions from the lysogenic to anti-immune states could be extracted by tracing back trial runs 
from A^ai back to Niyg. An analysis of these trajectories could be used to investigate the transition 
mechanism. However, in this work we have confined ourselves to computing the rate constant k. 

The frequency of transitions from the anti-immune to lysogenic states can be computed in the 
same way - here, our order parameter is simply taken to be 

A^' = total number of CI — total number of Cro (9) 
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To obtain the results presented in this paper, we used 8-85 interfaces, with 1000-10000 config- 
urations collected at A'^o and 500-10000 trial runs (precise values varied in each calculation). In 
addition, we averaged our results over 10-100 FFS calculations, to obtain an estimate for the error 
bars on our calculated rate constants. The final rate constants obtained are not dependent on the 
exact values of these parameters 



23l ? ]. For implementations of FFS with very few (~ 10 — 50) con- 
figurations at each interface, sampling errors may arise if the configurations most likely to progress 
to the final state are missed, especially if the order parameter poorly represents the progress of the 
transition 
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64] . These can be detected by the appearance of large variations between the rate 
constants calculated in different FFS runs. However, we have found this not to be the case when 
the chosen order parameter is reliable and the number of configurations is more than about 100, 



as in this work. 



X. SIMULATION RESULTS 



In this section, we present data on steady-state CI and Cro concentrations in the lysogenic and 
anti-immune states (IX Ap . depletion of free CI by nonspecific DNA binding, promoter activity as 
a function of lysogenic CI concentration (|X C\\ . switch stability for the Little mutants, calculated 
using FSS (jXPp and the effects of increasing the RNAp concentration (jXEp . 



A. Steady-state protein concentrations 

Figure [5] shows the average number of CI and Cro molecules in our simulations of the lysogenic 
and anti-immune states respectively, for different values of the parameters AGioop and AGnsb- 
The model produces steady-state protein concentrations in agreement with measured values. The 
number of CI molecules in a wild-type lysogen has been measured to be 220 [41[, but is highly 



variable uM- The concentration of Cro in the anti- immune state has been measured as ~ 400nM 
[1^ . so that in our assumed volume of 2^m^ we expect ~ 500 molecules. Our model predicts 
that the steady-state CI level in the lysogen is rather stable to changes in the nonspecific DNA 
binding strength AGnsb- However, if the looping interaction is absent, the lysogenic state cannot 
be sustained when AGnsb becomes large. 

Increasing the strength of the DNA looping interaction decreases the lysogenic CI concentration 



by a factor of up to 2, in agreement with the observations of Dodd et al 



8|. In contrast, the level 



of Cro in the anti-immune state increases strongly with the strength of the nonspecific binding 
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interaction, although it is insensitive to the DNA looping strength (as the DNA loop does not form 
in the anti-immune state). Pr is repressed by Cro binding at OrI or Or2. Depletion of free Cro 
by nonspecific binding reduces this autorepression, increasing expression of Pr, and increasing the 
total cellular Cro level. 



B. Depletion of free transcription factors by nonspecific DNA binding 



AGnsb ( kcal/mol) 


free monomers 


free dimcrs 


NSB dimers 


percentage free CI 


-2.1 


76 


75 


37 


75% 


-2.8 


52 


48 


76 


49% 


-3.5 


23 


24 


115 


24% 


-4.1 


18 


9 


132 


12% 



TABLE II: Partioning of the intracellular CI pool into free monomers, free dimers and nonspecifically bound 
dimers, for several values of AGnsb, assuming a total of 300 CI molecules in the cell and ignoring specific 
DNA binding reactions. The percentage of the CI pool which is free [not DNA-bound] is shown in the 
right-hand column. 

Table HI] shows the partioning of the total CI in the cell into free monomers, free dimers and 
nonspecifically bound dimers, assuming a total of 300 molecules in the cell and ignoring specific 
DNA binding. This demonstrates that the amount of free CI is dramatically depleted by nonspecific 
DNA binding, for the range of interaction strengths considered in this work. 



C. Promoter Activity 

In recent experiments by Dodd et al [3, Q], the activity of the Prm and Pr promoters was 
measured as a function of the intracellular CI concentration (supplied from a plasmid), for con- 

P Q 

structs with and without the DNA looping interaction [7|, [0]. Fitting the resulting data with a 
thermodynamic binding model using in vitro measured binding constants showed that CI binding 
to Or was much weaker than expected [3, Similar effects are well-known for other promoters, 



including lac 



42l |. and have been attributed to nonspecific binding of transcription factors to ge- 



nomic DNA. Dodd et al 



and other authors 



14 



631 ] fitted the measured promoter activity 



data with a model including nonspecific DNA binding. Although these promoter activity curves 
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may be distorted due to plasmid copy number variation [? ], the reduction in the apparent free CI 
concentration, apparently due to nonspecific binding, appears to be real and significant. 

400 r 




10000 



total 



FIG. 6: Computed activity of the Prm promoter, as a function of total cellular CI concentration. The red 
lines are computed with AGnsb = — 4.1kcal/mol; the blue lines without nonspecific DNA binding. The 
symbols show the data of Dodd et al, 2001 (circles) and Dodd et al, 2004 (squares), (a): Pum activity in 
absence of DNA looping interaction. 




10000 



total 



FIG. 7: Computed activity of the Prm promoter, as a function of total cellular CI concentration. The red 
lines are computed with AGnsb = — 4.1kcal/mol; the blue lines without nonspecific DNA binding. The 
symbols show the data of Dodd et al, 2001 (circles) and Dodd et al, 2004 (squares). Pbm activity, with 
DNA looping. The red and blue lines are computed with AGioop = — Ikcal/mol. The green line is computed 
with AGioop = — 3.7kcal/mol . The simulation data was multiphed by a scahng factor of 10000. 



In Figures [6] and [TJ we measure Prm activity as a function of fixed total CI concentration in 
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our model (CI levels are fixed by turning off protein production; no Cro is present). Increasing the 
nonspecific binding strength (magnitude of AGnsb) shifts the promoter activity curves towards 
ligher CI concentrations. The promoter activity data measured by Dodd et al, extracted from Refs 
3] and Q], are also shown in Figures [6] and [71 To convert Wild-type Lysogenic Units (WLU) to CI 
concentrations, we assumed a lysogenic CI concentration of 370nM [8|. Without nonspecific DNA 
binding, our model completely fails to agree with the measured Prm activity curves (blue lines). 
By fitting their data to a model with strong nonspecific DNA binding [AGnsb = — 6.2kcal/mol, 
with 6.76 X 10~'^M nonspecific binding site concentration], Dodd et al obtained a best-fit value 
of — 0.5kcal/mol for the DNA looping strength. In our model, the nonspecific binding is weaker 
compared to Dodd et al (AGnsb ~ 2 — 4kcal/mol; 1.0 x 10^ sites in a volume of 2;um - in the same 
volume, Dodd et al would have 1.35 x 10'' sites), and the DNA looping interaction is stronger. 
Dodd et al did not model the dynamics of the system, and therefore were not aiming to reproduce 
the stability of the lysogenic state. Our model, by contrast, is dynamical. We find that a strong 
looping interaction is needed to explain the stability and robustness of the lysogenic state; however, 
with our parameters, we do not fit the Dodd et al data well, as shown, for example, by the green 
line in Figure [71 Plasmid copy number variation may account for some of the discrepancy, as 
may our smaller value of AGnsb- However, the apparent discrepancy between the large value of 
^G^ioop needed to obtain agreement with experimental observations in dynamical simulations, and 
the smaller value obtained from fitting promoter activity data, remains to be resolved. 



D. Spontaneous switching times for the Little mutants 

We have used Forward Flux Sampling to compute spontaneous switching rates from the lyso- 
genic to the anti-immune state, and vice versa, for the Little mutants Or{121) and Or(323). The 
computed rates were converted into typical numbers of cell generations for which we expect these 
states to be stable, using a generation time of 34 min. The results are shown in IIIIl for a range of 
values of the parameters AGioop and AGnsb- 

We note that in modeling these mutants, one must decide which pair of CI dimers bound at Or 
experiences a cooperative binding interaction. In our simulations, we assumed cooperative binding 
interactions for the first two adjacent CI dimers to bind. 
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Mutant 


AGloop 


AGnsb 


Switching time 


Switching time 




kcal/mol 


kcal / mol 


lys. anti-imm. 


anti-imm. lys. 








(generations) 


(generations) 


Ofl(121) 


-5.2 


-4.1 


(1.1 ±0.4) X 10^° 


9.5 ±0.1 


Or{121) 


-3.7 


-2.8 


(4 ± 3) X 1015 


6.12 ±0.02 


Or{121) 


-1.0 


no NSB 


lysogen monostable 


anti-imm. unstable 


Oi? (323) 


-5.2 


-4.1 


450 ± 40 


450 ± 30 


Oij(323) 


-3.7 


-2.8 


(3 ±2) X lO^-i 


190 ± 20 


Ojj(323) 


-1.0 


no NSB 


130 ± 10 


(2.5 ±0.6) X 10^° 



TABLE III: Spontaneous switching times (inverse of calculated switching rates) for the Little mutants 
Oij(121) and Oij(323), computed as in Table 1 of the main text and assuming a cell generation time of 34 
min. 

E. Effects of RNA polymerase concentration 



We tested the effect of the free RNA polymerase concentration, which was fixed at 50nM in 
our model. This value was obtained by McClure, by fitting in vivo promoter expression data 
but has recently been challenged 



64l ]. We varied the free RNAp concentration and measured the 
total number of CI molecules in the lysogenic state. The results are plotted in Figure El Changing 
the free RNAp concentration can change the lysogenic CI concentration by about a factor of 
two. Although this parameter should be investigated further in future work, its value (within a 
reasonable physiological range) does not affect the qualitative features of our results. In fact, a 
higher free RNAp concentration is likely to further increase the stability of the lysogenic state. 



XI. SENSITIVITY OF OUR CONCLUSIONS TO PARAMETER CHOICE 

Our computational model involves several hundred chemical reactions. Many of the rate con- 
stants for these reactions are identical, so that the number of parameters is significantly less than 
the number of reactions. Nevertheless, parameter choice remains an important issue. We have 
taken most of the parameters for our model from the extensive biochemical data on the phage A 
switch available in the literature, as discussed in the main text. However, these values are subject 
to uncertainty, and in a few cases were unavailable, forcing us to make educated guesses (as detailed 
in the main text). It is therefore appropriate to investigate how sensitive our conclusions are to the 
choice of parameters. It is clearly impossible for us to prove that our conclusions hold throughout 
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FIG. 8: Effect of changing the free RNA polymerase concentration. The total number of CI molecules in 
the lysogenic state is plotted as a function of the free RNAp concentration. 



the many-dimensional space of all possible parameter variations. Nevertheless, we can strengthen 
our conclusions by testing the effects of several parameter changes, selected with the biology of the 
system in mind. To this end, we chose two parameters which are particularly uncertain, and are 
especially likely to influence the switching rate. We repeated our simulations for different values of 
these parameters. These are (1) the number of free RNA polymerase molecules in the cell, and (2) 
the average number of Cro molecules produced per mRNA molecule (translational burst size). The 
free RNAp concentration has not been measured directly in experiments, and inferred values are 
subject to a large range of uncertainty. The translational burst size is also rather uncertain from 
biochemical data. The burst size for Cro is likely to be particularly important because is possible 
that a burst of Cro molecules produced from a single mRNA might be sufficient to flip the switch. 
We varied the burst size by increasing the translation rate while proportionately decreasing the 
transcription rate, so as to keep the average Cro level in the cell approximately constant. Figures 
[9l and [TOl show the range of bistability for the model switch, as a function of the nonspecific DNA 
binding and DNA looping strengths, for free RNAp concentration doubled compared to the "stan- 
dard" model [lOOnM instead of 50nM] [9]and Cro translational burst size doubled [40 per transcript 
instead of 20] [TOl Both these plots show the same qualitative behaviour as we obtained for the 
"standard" parameter set, as shown in Fig. [3] of the main text. As the DNA looping strength 
increases, the range of nonspecific binding strengths over which the switch is bistable increases. 
We also used Forward Flux Sampling to compute switching rates for the lysogenic to anti-immune 
transition, in the case where the translational burst size for Cro was doubled. The results are 
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shown in Table IIVI 



no loop 
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FIG. 9: Range of bistability of the model switch as a function of the nonspecific DNA binding and DNA 
looping strengths, for a parameter set where the free RNAp concentration has been doubled [lOOnM instead of 
50nM] . Comparing these results to those of Fig. 3 in the main text shows that our main conclusion, namely 
that DNA looping stabilises the switch against nonspecific DNA binding, also holds for these modified 
parameter sets. 



Our main conclusions, that the lysogen is destabilized by nonspecific DNA binding but stabilized 
by the DNA looping interaction, remain the same with this modified parameter set. Comparing the 
two results without DNA looping in Table ITVl the switching time in the presence of nonspecific DNA 
binding is two orders of magnitude faster than when nonspecific binding is absent. Comparing the 
results with and without DNA looping, we find that the looping interaction stabilises the lysogen 
by a factor of 10^ in the absence of nonspecific binding and by a factor of 10^ when nonspecific 
binding is present. These switching rates, together with the bistability analysis of Figs. [9] and [101 
suggest that the conclusions arising from this work are likely to be robust to uncertainties in the 
parameter values. 
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FIG. 10: Range of bistability of the model switch as a function of the nonspecific DNA binding and DNA 
looping strengths, for a parameter set where the average number of Cro molecules produced per mRNA 
transcript has been doubled [40 instead of 20] by doubling the translation rate while halving the transcription 
rate. Comparing these results to those of Fig. [3] in the main text shows that our main conclusion, namely 
that DNA looping stabilises the switch against nonspecific DNA binding, also holds for these modified 
parameter sets. 



XII. COMPONENTS AND REACTIONS 



A. List of Components 

This is the list of components for our basic model, without non-specific binding and without 
looping. With the notation O(XYZ) we refer to the operator Or where X is bound to the site 
Or3, Y is bound to Or,2 and Z is bound to Ori. {X,Y,Z} = {0, R, C,Rp}, where corresponds 
to an empty site, R corresponds to CI dimer bound, C to a Cro dimer bound, and Rp to an RNA 
polymerase molecule bound. The Cro promoter, Pr overlaps with Or3, whereas the CI promoter, 
Prai overlaps with both Or2 and Ori, in our notation, RNA polymerase binds either to the X slot 
or to Y and Z simultaneously. 

1 CI 

2 Cro 
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AGioop 


AGnsb 


Switching time 


kcal/mol 


kcal/mol 


lys. — > anti-imm. 






(generations) 


no loop 


no NSB 


(1.5 ±0.2) X 10^ 


no loop 


-2.1 


(2.3 ±0.1) X 10^ 


-2.4 


no NSB 


(1.5 ±0.4) X 10" 


-2.4 


-2.1 


(2.7 ±0.2) X 1013 



TABLE IV: Spontaneous switching times (inverse of calculated switching rates) for the model switch with 
Cro translational burst size doubled, by doubling the translation rate while halving the transcription rate. 

3 CI2 

4 Cro2 

5 Rp 

6 O(OOO) 

7 O(ROO) 

8 O(ORO) 

9 O(OOR) 

10 O(COO) 

11 O(OCO) 

12 O(OOC) 

13 O(RpOO) 

14 O(ORp) 

15 O(RRO) 

16 O(ROR) 

17 O(RCO) 

18 O(ROC) 

19 O(RRp) 

20 O(ORR) 

21 O(CRO) 

22 O(ORC) 

23 0(RpRO 

24 O(COR) 

25 O(OCR) 
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26 O(RpOR) 

27 O(CCO) 

28 O(COC) 

29 O(CRp) 

30 O(OCC) 

31 O(RpCO) 

32 O(RpOC) 

33 O(RRR) 

34 O(RRC) 

35 O(RCR) 

36 O(CRR) 

37 O(RpRR) 

38 O(RCC) 

39 O(CRC) 

40 O(CCR) 

41 O(RpCR) 

42 O(RpRC) 

43 O(CCC) 

44 O(RpCC) 

45 O(RpRp) 

46 MCI 

47 MCro 

Adding nonspecific DNA binding requires 3 extra species: 

48 D 

49 DCI2 

50 DCro2 

When DNA looping is to be included in the model, we consider also the left operator Ol, which 
also has three binding sites for transcription factors. RNA polymerase can bind to promoter Pl, 
which overlaps with the binding site Oli- The two operator can interact via a DNA loop stabilised 
by the octamerization of two CI tetramers which are already bound on adjacent sites on the two 
operators. We label these looped states according to which pairs of CI dimers form the octamer. 
Looped states obtained by interaction between CI tetramers bound to Ori, Or2jOli and Ol2 
are denoted OLRl. States 0LR2 have loops formed by interactions between tetramers bound 
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to Ori, Or2 5 Ol2 and Ols- States denoted 0LR3 have loops formed by interaction between 
tetramers bound to 0^2, Or3, Oli and Ol2- Finally, states denoted 0LR4 have loops formed by 
interaction between tetramers bound to Or2, Or3, Ol2 and Ols- Each of these loop categories 
has two additional sites to which transcription factor dimers can bind. We therefore denote a 
particular looped configuration by OLRX(YZ), where X refers to loop type 1-4 (as above), Y 
denotes which species occupies the non-octamerized site on Ol and Z labels the occupation of the 
non-octamerized site on Or. The extra species required in the model with looping are: 

51 OL(OOO) 

52 OL(ROO) 

53 OL(ORO) 

54 OL(OOR) 

55 OL(COO) 

56 OL(OCO) 

57 OL(OOC) 

58 OL(OORp) 

59 OL(RRO) 

60 OL(ROR) 

61 OL(RCO) 

62 OL(ROC) 

63 OL(RORp) 

64 OL(ORR) 

65 OL(CRO) 

66 OL(ORC) 

67 OL(ORRp) 

68 OL(COR) 

69 OL(OCR) 

70 OL(CCO) 

71 OL(COC) 

72 OL(CORp) 

73 OL(OCC) 

74 OL(OCRp) 

75 OL(RRR) 

76 OL(RRC) 



77 OL(RRRp) 

78 OL(RCR) 

79 OL(CRR) 

80 OL(RCC) 

81 OL(RCRp) 

82 OL(CRC) 

83 OL(CRRp) 

84 OL(CCR) 

85 OL(CCC) 

86 OL(CCRp) 

87 OLRl(OO) 

88 OLRl(OR) 

89 OLRl(OC) 

90 OLRl(ORp) 

91 OLRl(RO) 

92 OLRl(RR) 

93 OLRl(RC) 

94 OLRl(RRp) 

95 OLRl(CO) 

96 OLRl(CR) 

97 OLRl(CC) 

98 OLRl(CRp) 

99 OLR2(00) 

100 OLR2(0R) 

101 OLR2(0C) 

102 OLR2(0Rp) 

103 OLR2(R0) 

104 0LR2(RR) 

105 0LR2(RC) 

106 0LR2(RRp) 

107 OLR2(C0) 

108 0LR2(CR) 

109 0LR2(CC) 



110 0LR2(CRp) 

111 OLR3(00) 

112 OLR3(0R) 

113 OLR3(0C) 

114 OLR3(R0) 

115 0LR3(RR) 

116 0LR3(RC) 

117 OLR3(C0) 

118 0LR3(CR) 

119 0LR3(CC) 

120 OLR4(00) 

121 OLR4(0R) 

122 OLR4(0C) 

123 OLR4(R0) 

124 0LR4(RR) 

125 0LR4(RC) 

126 OLR4(C0) 

127 0LR4(CR) 

128 0LR4(CC) 

129 OLR2(RpO) 

130 0LR2(RpR) 

131 0LR2(RpC) 

132 0LR2(RpRp) 

133 OLR4(RpO) 

134 0LR4(RpR) 

135 0LR4(RpC) 



B. List of reactions 

This set is obtained with AGioop = —5.2 kcal/mol and AGnsb = —3.5 kcal/mol. 

1) CI + CI ^ CI2 k = 0.628319 

2) CI2 ^ CI + CI k = 5.705413 



3) Cro + Cro Croa k = 0.628319 

4) Cro2 ^ Cro + Cro k = 280.199086 

5) O(OOO) + CI2 ^ O(ROO) k = 0.314159 

6) 0(ROO ^ O(OOO) + CI2 k = 38.256936 

7) O(OOO) + CI2 ^ O(ORO) k = 0.314159 



8) 


O(ORO) ^ O(OOO) + CI2 k = 7.551827 


9) O(OOO) + CI2 ^ O(OOR) k = 0.314159 


10) 


O(OOR) O(OOO) + CI2 k = 


0.294263 


11) 


O(OOO) + Cro2 O(COO) k = 


= 0.314159 


12) 


O(COO) O(OOO) + Cro2 k = 


= 0.068319 


13) 


O(OOO) + Cro2 O(OCO) k = 


= 0.314159 


14) 


O(OCO) O(OOO) + Cro2 k = 


= 4.641460 


15) 


O(OOO) + Cro2 O(OOC) k = 


= 0.314159 


16) 


O(OOC) O(OOO) + Cro2 k = 


= 0.662315 


17) 


O(OOO) + Rp ^ O(RpOO) k = 


= 0.314159 


18) 


O(RpOO) O(OOO) + Rp k = 


= 1.490712 


19) 


O(OOO) + Rp ^ O(ORp) k = 


0.314159 


20) 


O(ORp) O(OOO) + Rp k = 


0.294263 


21) 


O(ROO) + CI2 O(RRO) k = 


= 0.314159 


22) 


O(RRO) O(ROO) + CI2 k = 


= 0.068319 


23) 


O(ROO) + CI2 O(ROR) k 


= 0.314159 


24) 


O(ROR) O(ROO) + CI2 k 


= 0.294263 




Ul^xlUUJ + '^r02 — > U(^rlLyUj K 


n Q 1 /1 1 KQ 

— u.oi4ioy 


26) 


O(RCO) ^ O(ROO) + Cro2 k 


= 4.641460 


27) 


O(ROO) + Cro2 ^ O(ROC) k 


= 0.314159 


28) 


O(ROC) ^ O(ROO) + Cro2 k 


= 0.662315 


29) 


O(ROO) + Rp ^ O(RRp) k = 


= 0.314159 


30) 


O(RRp) ^ O(ROO) + Rp k = 


= 0.294263 


31) 


O(ORO) + CI2 ^ O(RRO) k = 


= 0.314159 


32) 


O(RRO) ^ O(ORO) + CI2 k = 


= 0.346100 


33) 


O(ORO) + CI2 ^ O(ORR) k = 


= 0.314159 


34) 


O(ORR) ^ O(ORO) + CI2 k = 


= 0.003683 


35) 


O(ORO) + Cro2 ^ O(CRO) k 


= 0.314159 
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;rco) + 


CI2 ^ O(RCR) k = 


n Qi /1 1 c^n 
U.oi4ioy 


iizj 


U 


;rcr) - 


> O(RCO) + CI2 k = 


u.zy4zoo 


iioj 


U 


^RCO) + 


Cro2 O(RCC) k = 


= U.oi4ioy 


ii4J 


U 


;rcc) - 


> O(RCO) + Cro2 k = 


= U.ioU t oy 


iioJ 


U 


;roc) + 


CI2 O(RRC) k = 


n Q 1 /1 1 
U.oi4ioy 


ilDj 


r\ 
U 


^RRC) - 


> O(ROC) + CI2 k = 


U.UOooiy 


1 1 7^ 
ii ( J 


(J 


^ROC) + 


Cro2 O(RCC) k = 


— n QI /1 1 c^o 

- u.oi4ioy 


iioJ 


(J 


;rcc) - 


> O(ROC) + Cro2 k = 


1 vr^QQl A 
- i./i300i4 


iiyj 


U 


^CRO) + 


CI2 O(CRR) k = 


U.oi4ioy 


1 o^^ 
i2U J 


(J 


^CRR) - 


> O(CRO) + CI2 k = 


U.UUoDoo 


LZL j 


(J 


^CRO) + 


Cro2 O(CRC) k = 


= U.oi4ioy 


izz J 


U 


;cRC) - 


> O(CRO) + Cro2 k = 


= U.DDzoiO 


i2c) J 


(J 


^COR) + 


CI2 O(CRR) k = 


n Q 1 /1 1 c^n 
U.oi4ioy 


iz4 J 




;CRR) - 


> O(COR) + CI2 k = 


u.uy4ouy 


i2o J 


u 


;COR) + 


Cro2 O(CCR) k = 


= U.oi4ioy 


1 o^^^ 
i2D J 


(J 


;ccR) - 


> O(COR) + Cro2 k = 


1 VF^QQI /I 

= i. / 0ooi4 




u 


'OCR) + 


CI2 O(RCR) k = 


n Q 1 /1 1 
U.oi4ioy 


128 J 


(J 


;rcr) - 


> O(OCR) + CI2 k = 


oo.zooyoD 


1 o^^ 

i2y J 


(J 






n QI /1 1 F^n 
= U.oi4ioy 


ioU J 


p\ 
(J 


:ccR) - 


> O(OCR) + Cro2 k = 


= U.UzooUo 


131) 





;ocr) + 


Rp ^ O(RpCR) k = 


= 0.314159 


132) 





;RpCR) 


O(OCR) + Rp k = 


= 1.490712 


133) 





;ORC) + 


CI2 ^ O(RRC) k = 


0.314159 


134) 





;rrc) - 


> O(ORC) + CI2 k = 


0.346100 



135) O(ORC) + Cro2 ^ O(CRC) k = 0.314159 

136) O(CRC) ^ O(ORC) + Croa k = 0.068319 

137) O(ORC) + Rp ^ O(RpRC) k = 0.314159 

138) O(RpRC) ^ O(ORC) + Rp k = 1.490712 

139) O(CCO) + CI2 ^ O(CCR) k = 0.314159 

140) O(CCR) ^ O(CCO) + CI2 k = 0.294263 

141) O(CCO) + Cro2 ^ O(CCC) k = 0.314159 

142) O(CCC) ^ O(CCO) + Cro2 k = 0.407068 

143) O(COC) + CI2 ^ O(CRC) k = 0.314159 

144) O(CRC) ^ O(COC) + CI2 k = 7.551827 

145) O(COC) + Cro2 ^ O(CCC) k = 0.314159 

146) O(CCC) ^ O(COC) + Cro2 k = 1.077612 

147) O(OCC) + CI2 ^ O(RCC) k = 0.314159 

148) O(RCC) ^ O(OCC) + CI2 k = 38.256936 

149) O(OCC) + Cro2 ^ O(CCC) k = 0.314159 

150) O(CCC) ^ O(OCC) + Cro2 k = 0.080354 

151) O(OCC) + Rp ^ O(RpCC) k = 0.314159 

152) O(RpCC) ^ O(OCC) + Rp k = 1.490712 

153) O(RpRO) + CI2 ^ O(RpRR) k = 0.314159 

154) O(RpRR) ^ O(RpRO) + CI2 k = 0.003683 

155) O(RpRO) + Cro2 ^ O(RpRC) k = 0.314159 

156) O(RpRC) ^ O(RpRO) + Cro2 k = 0.662315 

157) O(RpOR) + CI2 ^ O(RpRR) k = 0.314159 

158) O(RpRR) ^ O(RpOR) + CI2 k = 0.094509 

159) O(RpOR) + Cro2 ^ O(RpCR) k = 0.314159 

160) O(RpCR) ^ O(RpOR) + Cro2 k = 4.641460 

161) O(RpCO) + CI2 ^ O(RpCR) k = 0.314159 

162) O(RpCR) ^ O(RpCO) + CI2 k = 0.294263 

163) O(RpCO) + Cro2 ^ O(RpCC) k = 0.314159 

164) O(RpCC) ^ O(RpCO) + Croa k = 0.130739 

165) O(RpOC) + CI2 ^ O(RpRC) k = 0.314159 

166) O(RpRC) ^ O(RpOC) + CI2 k = 7.551827 

167) O(RpOC) + Cro2 ^ O(RpCC) k = 0.314159 
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iuo J 


O(RpCC) ^ O(RpOC) + Cro2 k = 1.753314 


ioy J 


O(RpOO) 


O(OOO) + Rp + MCI k = 


0.001000 


i /UJ 


O(RpRO) - 


> O(ORO) + Rp + MCI k = 


= 0.011000 


1 71 ^ 


O(RpOR) - 


> O(OOR) + Rp + MCI k = 


= 0.001000 


1 

i (Zj 


O(RpCO) - 


> O(OCO) + Rp + MCI k = 


= 0.001000 


1 7Q^ 
i /oj 


O(RpOC) - 


> O(OOC) + Rp + MCI k = 


= 0.001000 




O(RpRR) - 


O(ORR) + Rp + MCI k 


= 0.011000 




O(RpRC) - 


O(ORC) + Rp + MCI k 


= 0.011000 


1 7^^ 

i /Dj 


O(RpCR) - 


-y O(OCR) + Rp + MCI k 


= 0.001000 


1 77^ 
III) 


O(RpCC) - 


O(OCC) + Rp + MCI k 


= 0.001000 


1 7Q^ 

i (Oj 


O(RpRp) - 


^ O(ORp) + Rp + MCI k 


= 0.001000 


1 7^^ 


O(ORp) ^ 


O(000)+ Rp + MCro k = 


014000 


1 o^\^ 

loUj 


O(RRp) ^ 


O(ROO) + Rp + MCro k = 


= 014000 


lolj 


O(CRp) ^ 


O(COO) + Rp + MCro k = 


= 014000 

• KJ -M- KJ KJ KJ 


18zj 


O(RpRp) - 


■> O(RpOO) + Rp + MCro 


k = 014000 

J.V KJ * KJ -L- \J KJ \J 


1 oo\ 

io6) 


MCI ^ MCI + CI k = 0.034656 




184j 


MCro MCro + Cro k = 0.115520 




loo) 


MCI ^ k 


= 0.005776 




Im) 


MCro ^ k = 0.005776 




187) 


CI ^ k = 


0.000340 




188) 


Cro ^ k 


= 0.000606 




189) 


CI2 ^ k : 


= 0.000340 




190) 


Cro2 ^ k 


= 0.000340 





If the system includes looping, the reaction schemes should be extended with the following reactions 
(rate constants calculated for AGioop = —5.2 kcal/mol): 

191) OL(OOO) + CI2 OL(ROO) k = 0.314159 

192) OL(ROO) OL(OOO) + CI2 k = 0.346100 

193) OL(OOO) + CI2 ^ OL(ORO) k = 0.314159 

194) OL(ORO) ^ OL(OOO) + CI2 k = 0.563117 

195) OL(OOO) + CI2 ^ OL(OOR) k = 0.314159 

196) OL(OOR) ^ OL(OOO) + CI2 k = 0.035701 

197) OL(OOO) + Cro2 OL(COO) k = 0.314159 
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1 OQ\ 

lyo J 


OL(COO) ^ 


OL(OOO) + Cro2 k = 


U.UDooiy 


lyy ) 


OL(OOO) + Cro2 ^ OL(OCO) k = 


U.oi4ioy 


zUU J 


OL(OCO) ^ 


OL(OOO) + Cro2 k = 


4.D4i4DU 


zUi J 


OL(OOO) + Cro2 ^ OL(OOC) k = 


n Q 1 /1 1 
U.oi4ioy 


zUz J 


OL(OOC) 


OL(OOO) + Cro2 k = 


U.DOZOiO 


zUo J 


OL(OOO) + 


Rp OL(OORp) k = 


n Q 1 /1 1 c^n 
U.oi4ioy 


OH/I ^ 
zU4 J 


OL(OORp) - 


OL(OOO) + Rp k = 


Z.UbZi (0 


zUo J 


OL(ROO) + 


CI2 OL(RRO) k = 


n Q 1 /1 1 
U.oi4ioy 


zUo J 


OL(RRO) - 


> OL(ROO) + CI2 k = 


u.uuy r4y 


zU /^J 


OL(ROO) + 


CI2 OL(ROR) k = 


n Q 1 /1 1 
U.oi4ioy 


zUo J 


OL(ROR) - 


^ OL(ROO) + CI2 k = 


U.UoO / Ui 


zuy J 


OL(ROO) + 


Cro2 OL(RCO) k 


= U.oi4ioy 


ziU J 


OL(RCO) - 


OL(ROO) + Cro2 k 


= 4.D4i4oU 


ziij 


OL(ROO) + 


Cro2 OL(ROC) k = 


= U.oi4ioy 


ziz J 


OL(ROC) - 


> OL(ROO) + Cro2 k = 


= U.DDzoit) 


Zio) 


OL(ROO) + 


Kp OL(RORp) k = 


- n Qi /1 1 


OI /I ^ 
Zl^j 


OL(RORp) 


OL(ROO) + Rp k = 


- z.Uozi 1 


OI 

zio J 


OL(ORO) + 


CI2 OL(RR0 k = 


U.oi4ioy 


OI 

zib J 


OL(RRO) - 


* OL(ORO) + CI2 k = 


u.uuoyyz 


OI '7^ 

Zi (^J 


OL(ORO) + 


CI2 OL(ORR) k = 


n Q 1 /1 1 
U.oi4ioy 


OI Q^ 
zioj 


OL(ORR) - 


> OL(ORO) + CI2 k = 


U.UUUOio 


OI 

ziy J 


OL(ORO) + 


Cro2 OL(CRO) k = 


n Q 1 /1 1 c^n 
= U.oi4ioy 


zzU J 


OL(CRO) - 


^ OL(ORO) + Cro2 k = 


= U.Ubooiy 


001 ^ 
zzl J 


OL(ORO) + 


Cro2 OL(ORC) k = 


= U.oi4ioy 


ooo^ 
zzz J 


OL(ORC) - 


> OL(ORO) + Cro2 k = 


= U.DDZOiO 


ZZO J 


OT /HRn^ _l_ 
(Jij[Ut\X}) -\- 


xip — > Wij(^Urlxipj K - 


= u.oi4ioy 


00/1 ^ 
zz4 J 


OL(ORRp) 


OL(ORO) + Rp k = 


= Z.UOZi (0 


ZZO J 


OL(OOR) + 


CI2 ^ OL(ROR) k = 


n Q 1 /1 1 


ZZD J 


OL(ROR) - 


> OL(OOR) + CI2 k = 


n Q/i fii nn 


227) 


OL(OOR) + 


CI2 ^ OL(ORR) k = 


0.314159 


228) 


OL(ORR) - 


* OL(OOR) + CI2 k = 


0.009749 


229) 


OL(OOR) + 


Cro2 ^ OL(COR) k = 


= 0.314159 


230) 


OL(COR) - 


> OL(OOR) + Cro2 k = 


= 0.068319 
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OQ1 \ 

zoi J 


Ulj^UUxlj + 


L/ro2 — > U-HUL/Kj K = 


n Q 1 /1 1 F:n 
= U.oi4iOy 


zoz J 


UL{()kjix) — 




A P.A 1 AP.r\ 

- 4.t)4i4t)U 


zoo ) 






U.oi4iOy 


zo4 J 


Uh[(^tW) — 


^ p>T ('P'nnA 1 PT 1^ 


U.OOoii t 


zoo ) 




PT ^ p>T (r^r\T} \ ^r 


n Q 1 /1 1 c^n 
U.oi4iOy 


zoo J 


UL[Ly\j\x) — 


^ p>T fr^r\r\\ 1 pt 1^ 


U.UoO (Ui 


ZO ( J 




Pt.^ ^ p>T /'P'P'n\ 1, 
L/r02 — > *Jlj(L/L/Uj K = 


n Q 1 /1 1 
= U.oi4ioy 


zoo J 


UL[LjLj\j) — 


> Ulj(L/UUj + L/r02 K = 


- i. /Oooi4 


zoy J 




Pt.^ ^ p>T /'P'nP\ u 


n Q 1 /1 1 c^n 
= U.oi4ioy 


z4U J 


Uij(^UUL/j — 


K OT /'^~'nn^ 1 p».^ u 
> UL/(L/UUj + L'r02 K = 


= U.DDZOiO 


z4i J 


Uij(L/UUj + 


Kp — > Ulj(L/Urtpj K = 


: U.oi4iOy 


Z4z J 




V r^T fr^nri\ 1 t>>^ 1^ 
— > U1j(^L/UUJ + Kp K = 


: z.Ubzi rO 


z4o J 


Ulj(UL/Uj + 


PT V p>T /'^^P'o^ 1^ 


n Q 1 /1 1 
U.oi4ioy 


z44 J 




^ p>T /'nPriA 1 PT 1, 


n Q /I fii nn 
U.o4DiUU 


z4o J 




PT V p>T fnCD^ ^r 


U.oi4iOy 


Z4D J 




i OT lr\C*c\\ _L PT 1^ 

> Wlj(UL/Uj + Ui2 K — 


U.UoO ( Ui 


z4^J 




c^^r^ K p>T ('P'r~'n\ 1, 
L/r02 — > (Jlj(L/L/Uj K = 


n Q 1 /1 1 c^o 
= U.oi4ioy 


o/i Q^ 
z4o J 


*Jlj(L/L/Uj — 


> Ulj(UL/Uj + L/r02 K = 


= U.UzooUo 


O/I 

z4y J 




L/r02 — > (JL[U(^(^) K = 


n Q 1 /1 1 c^o 
= U.oi4ioy 


ZOU J 


*Jlj(^UL/Uj — 


» Ulj(ULyUj + L/ro2 K = 


= U.ioU/oy 


ZOi J 




Kp — !• UL/(UUKpj K = 


n Q 1 /1 1 c^n 
U.oi4iOy 


o p^o^ 
zoz J 


(JL[Uvjsxp ) 


— > (JLI^UUUJ + Kp K = 


: Z.Ubzi to 


zOo ) 


(Jlj^UUL/j + 


PT ^ P>T /Tsnp^ 1, 


n Q 1 /1 1 c^n 
U.oi4iOy 


z04 J 


Pit /'^^nP'^ 


^ p»T /'nnp'^ 1 pt 1^ 
> UL/I^UUL/j + Ui2 K = 


n Q /I 1 nn 
U.o4DiUU 


ZOO ) 


Uij(UUL/j + 


L/i2 — > Ulj(^UKL/j K = 


U.oi4iOy 


zOD ) 




> Ulj(UUL/j + L/I2 k; = 


U.Oboii / 


ZO ( J 




Pt.^ ^ PIT /'pnp\ u 


n Q 1 /1 1 c^n 
= U.oi4ioy 


zOo ) 


OL(COC) - 


> OL(OOC) + Cro2 k = 


n n^QQi n 
= U.UOooiy 


zoy J 


OL(OOC) + 


Cro2 ^ OL(OCC) k = 


n Q 1 /1 1 r^n 
= U.oi4ioy 


260) 


OL(OCC) - 


> OL(OOC) + Cro2 k = 


= 0.916213 


261) 


OL(OORp) + CI2 ^ OL(RORp) k 


= 0.314159 


262) 


OL(RORp) 


OL(OORp) + CI2 k 


= 0.346100 


263) 


OL(OORp) + CI2 ^ OL(ORRp) k 


= 0.314159 
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Z04J 


Ulj^UKxlpj — > Uij^UUxlpj + LjL2 k = 


= U.ODoii / 


ZOO ) 


Ulj^UUHpj + Ur02 — > (Jij^UUKpJ K 


= U.oi4ioy 


ZOO ) 


Uij^UUKpj — > Ulj^UUxlpj + U102 K 


= U.UDooiy 


ZD{ ) 


Ulj^UUxlpj + Ur02 — > Ulj^UUKpj K 


= U.oi4ioy 


zoo ) 


Ulj(^UUKpj — > Uh[U\jtxp) + Ur02 K 


A dA^ A fin 

= 4.D4i4DU 


zoy J 


Uljl^xlKUj + LjL2 — > UL[lxixtx) K = 


n Q 1 /1 1 c^n 
U.oi4ioy 


z / UJ 


Uljl^xtKxtj — > UL/I^KxtUj + L/i2 ii — 


U.Uoo ( Ui 


z / ij 


UL[txt\,U) + Ur02 — ^ U-UKKUj K = 


= U.oi4ioy 


Zi Z) 


Ulj^xlKUj — > UL[lxt\l)) + U102 K = 


- U.bDZOiO 


z /oj 


Ulj^xlKUj + xlp — > UL[lxixlxp) K = 


: U.oi4ioy 


z /4j 


Ulj^KKHpj — > UL[lxr\l)) + Kp K = 


z.Uozi /O 


z /Oj 


U-L^xlUxlj + LjL2 — > (JL[lxixlx) K = 


U.oi4ioy 


z /Dj 


U-L^xlKxlj — > Uij(KUKj + 012 K = 


u.uuy ( 4y 


Z / (^J 


Uljl^KUKj + Ur02 — > Ub^KUKJ K = 


= U.oi4ioy 


z / oj 


Ulj^xlLyxlj — > Uij^KUKj + '^102 K - 


- A dA 1 /ifin 

- 4.u4i4uU 


z / yj 


(JL(U]xix) + Lyi2 — > KJlj[ixsxlx) K — 


n Qi /1 1 c^o 
U.oi4ioy 


zoU J 


(Jlj^KKKj — > UHUKKJ + L/i2 K = 


U.o4DiUU 


zoi J 


U-HUKKj + L'r02 — > Uij(UKKj K = 


= U.oi4ioy 


zoz J 


(JIj^L/KKJ — > U1j(UKKJ + Ur02 K = 


= U.UDooiy 


zoo J 


(Jlj^KL/Uj + Ui2 — > Wlj^KUrtj K = 


U.oi4ioy 


Zo4 J 




U.UoO ( Ui 


zoo J 


(Jlj^KLyUj + Ur02 — > wL/l^xlUUj K : 


= U.oi4ioy 


zoo ) 


(Jlj^KLyUj — > UL[r{Xj(J) + Uro2 K : 


= U.ioUMy 


zo /^J 


Uij^ixUUj + Kp — > ULi^KUKpj K = 


= U.oi4ioy 


zoo ) 


(Jlj^xlL/Kpj — 5- Uij(^KLyUJ + Kp K = 


= Z.UbZi 10 


zoy ) 


(Jlj(KULyJ + Ui2 — > U1j(KKUJ k = 


U.oi4ioy 


zyu J 


Ulj^KKUj — > UL[txULj) + Ui2 K = 


U.uuy /4y 


zyi J 


*J1j(^KUL/J + Ur02 — ^ UL/^KULyj K = 


= U.oi4it)y 


zyz J 


(Jlj(KL/Uj — > UL[IxULj) + Lj1C02 = 


= i. 106614: 


293) 


OL(CRO) + CI2 ^ OL(CRR) k = 


0.314159 


294) 


OL(CRR) ^ OL(CRO) + CI2 k = 


0.000618 


295) 


OL(CRO) + Cro2 ^ OL(CRC) k = 


= 0.314159 


296) 


OL(CRC) ^ OL(CRO) + Cro2 k = 


= 0.662315 
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zy I J 


Uij^UilUj + 


Kp — > Uij^UKxlpj K = 


: U.oi4ioy 


zyo J 


Uij(L/ilKp j 


— > Uij(^UitUj + itp K — 


- z.UDzi (0 


zyy J 


Uij^UUKj + 


Ui2 — > Uij^UKKj K = 


U.oi4ioy 


oUU J 


Uij(L/ilKj — 


> UL/(UUilj + L/i2 K = 


U.UUy /4y 


oUi J 


Uij(^L/UKj + 


Ur02 — ^ Uij(UUKJ K = 


= U.oi4iOy 


oUz J 


Uij^L-L-Kj — 


> Uij(^UUKJ + U102 K = 


= i. ( 0ooi4 


oUo J 


Uij(^UUKj + 


Ui2 — > UL/(KUKJ K = 


U.oi4iOy 


oU4 J 


Uij^ilLyilj — 


> Uij^ULyilj + Ui2 K = 


U.o4DiUU 


aUb J 


Uij^ULyKj + 


U102 — > UL/(UUKj K = 


= U.oi4iOy 


oUD J 




> (JL[UK^ix) + UIO2 K - 


- U.UzOoUo 


6\}t ) 


Uij^UKUj + 


Ui2 — > UL[lxriLj) K = 


U.oi4iOy 


oUo J 


\J\^\ixt\Kj) — 


Uij(UilLvJ + Ui2 K = 


u.uuoyyz 


auy J 


Uij^UKUj + 


L'r02 — > UL[(^tiXj) K = 


= U.oi4iOy 


oiU J 


Uij^LyilUj — 


> Uij^UilUj + L/r02 K = 


= U.UOooiy 


oiij 


Uij^LvL/Uj + 


Ui2 — Uij^L-L/Kj K = 


U.oi4iOy 


OiZJ 


Oil tCOT) \ 


> kJL[(^Lj\J) + Ui2 K — 


U.Uoo ( Ui 


oio J 


Uij^LyUUj + 


Ur02 — > (Jij(L/UUJ K = 


= U.oi4iOy 


oi4 j 


UL[(^Lj(^ ) — 


> Uij(uuuj + 0102 K = 


= U.4U ( Udo 


oio J 


Uij(L/UUj + 


Kp — > Uij^UUKpj K = 


U.oi4ioy 


oiD J 


Uij(L/UKp ) 


— > Uij^UOUj + Kp K = 


z.Uozi (0 


Oi ( j 


Uij^UUUj + 


Ui2 — > \JIj[KjSXKj) K = 


U.oi4iOy 


oio J 




J' Uij(UUUJ + Ui2 K = 


U.Oboii / 


oiy J 


(Jij(L/UUj + 


r^v^ K c\i (r^oo\ 1, 


= u.oi4ioy 


ozU j 


Pit (oooy 


Uij(UUUj + 0102 K = 


= i.U ( (biz 


o2i J 


Uij^LH^Lyj + 


Ui2 — > Uij(KL/U J K = 


n Q 1 /1 1 F;n 
U.oi4iOy 


o2z J 


OL(RCC) - 


^ OL(OCC) + CI2 k = 


n Q /I ^ 1 nn 
U.o4DiUU 


ozo j 


OL(OCC) + 


Cro2 ^ OL(CCC) k = 


n Q 1 /1 1 F;n 
= U.oi4iOy 


oz4 J 


OL(CCC) - 


> OL(OCC) + Cro2 k = 


= 0.080354 


ozo J 


OL(ORRp) + CI2 ^ OL(RRRp) k 


= 0.314159 


326) 


OL(RRRp) 


OL(ORRp) + CI2 k 


= 0.005992 


327) 


OL(ORRp) + Cro2 ^ OL(CRRp) k = 0.314159 


328) 


OL(CRRp) 


OL(ORRp) + Cro2 k = 0.068319 


329) 


OL(RORp) + CI2 ^ OL(RRRp) k 


= 0.314159 



330) OL(RRRp) ^ OL(RORp) + CI2 k = 0.009749 

331) OL(RORp) + Cro2 ^ OL(RCRp) k = 0.314159 

332) OL(RCRp) ^ OL(RORp) + Cro2 k = 4.641460 

333) OL(OCRp) + CI2 ^ OL(RCRp) k = 0.314159 

334) OL(RCRp) ^ OL(OCRp) + CI2 k = 0.346100 

335) OL(OCRp) + Cro2 ^ OL(CCRp) k = 0.314159 

336) OL(CCRp) ^ OL(OCRp) + Cro2 k = 0.025808 

337) OL(CORp) + CI2 ^ OL(CRRp) k = 0.314159 

338) OL(CRRp) ^ OL(CORp) + CI2 k = 0.563117 

339) OL(CORp) + Cro2 ^ OL(CCRp) k = 0.314159 

340) OL(CCRp) ^ OL(CORp) + Cro2 k = 1.753314 

341) OL(ORR) + O(ORR) ^ OLRl(OO) k = 62.095095 

342) OLRl(OO) ^ OL(ORR) + O(ORR) k = 0.008046 

343) OL(ORR) + O(RRR) ^ OLRl(OR) k = 62.095095 

344) OLRl(OR) ^ OL(ORR) + O(RRR) k = 0.008046 

345) OL(ORR) + O(CRR) ^ OLRl(OC) k = 62.095095 

346) OLRl(OC) ^ OL(ORR) + O(CRR) k = 0.008046 

347) OL(ORR) + O(RpRR) ^ OLRl(ORp) k = 62.095095 

348) OLRl(ORp) ^ OL(ORR) + O(RpRR) k = 0.008046 

349) OL(RRR) + O(ORR) ^ OLRl(RO) k = 62.095095 

350) OLRl(RO) ^ OL(RRR) + O(ORR) k = 0.008046 

351) OL(RRR) + O(RRR) ^ OLRl(RR) k = 62.095095 

352) OLRl(RR) ^ OL(RRR) + O(RRR) k = 0.008046 

353) OL(RRR) + O(CRR) ^ OLRl(RC) k = 62.095095 

354) OLRl(RC) ^ OL(RRR) + O(CRR) k = 0.008046 

355) OL(RRR) + O(RpRR) ^ OLRl(RRp) k = 62.095095 

356) OLRl(RRp) ^ OL(RRR) + O(RpRR) k = 0.008046 

357) OL(CRR) + O(ORR) ^ OLRl(CO) k = 62.095095 

358) OLRl(CO) ^ OL(CRR) + O(ORR) k = 0.008046 

359) OL(CRR) + O(RRR) ^ OLRl(CR) k = 62.095095 

360) OLRl(CR) ^ OL(CRR) + O(RRR) k = 0.008046 

361) OL(CRR) + O(CRR) ^ OLRl(CC) k = 62.095095 

362) OLRl(CC) ^ OL(CRR) + O(CRR) k = 0.008046 
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ODO J 


OL(CRR) + O(RpRR) ^ OLRl(CRp) k = 62.095095 


/I ^ 
od4 j 


OLRl(CRp) ^ OL(CRR) + O(RpRR) k = 0.008046 


ODO ) 


OLRl(OO) + CI2 OLRl(RO) k = 


0.314159 


ODOj 


OLRl(RO) OLRl(OO) + CI2 k = 


0.346100 




OLR1(00)+ CI2 OLRl(OR) k = 


0.314159 


ODO J 


OLRl(OR) OLRl(OO) + CI2 k = 


38.256936 


ooy ) 


OLRl(OO) + Cro2 OLRl(CO) k = 


= 0.314159 


o / UJ 


OLRl(CO) OLRl(OO) + Cro2 k = 


= 0.068319 


Oil) 


OLRl(OO) + Cro9 OLRlfOC) k = 


= 314159 


O / Zj 


OLRl(OC) OLRl(OO) + Cro2 k = 


= 0.068319 




OLRl(OO) + Rp ^ OLRl(ORp) k = 


= 0.314159 




OLRl(ORp) OLRl(OO) + Rp k = 


= 1.490712 


o/oj 


OLRl(RO) + CI2 OLRl(RR) k = 


= 0.314159 


oiu) 


OLRl(RR) OLRl(RO) + CI2 k = 


= 0.294263 


O M J 


OLRl(RO) + Cro2 OLRl(RC) k 


= 0.314159 


olo) 


OLRl(RC) OLRl(RO) + Cro2 k 


= 0.068319 




OLRl(RO) + Rp ^ OLRl(RRp) k 


= 0.314159 


doU ) 


OLRl(RRp) OLRl(RO) + Rp k 


= 1.490712 




OLRl(CO) + CI2 OLRl(CR) k = 


= 0.314159 


ooz ) 


OLRl(CR) OLRl(CO) + CI2 k = 


= 38.256936 


ooo j 


OLRl(CO) + Cro2 OLRl(CC) k 


= 0.314159 


c)o4 J 


OLRl(CC) OLRl(CO) + Cro2 k 


= 0.068319 


ooo J 


OLRl(CO) + Rp ^ OLRl(CRp) k 


= 0.314159 


OOO ) 


OLRl(CRp) OLRl(CO) + Rp k 


= 1.490712 




OLRl(OR) + CI2 OLRl(RR) k = 


= 0.314159 


ooo ) 


OLRl(RR) OLRl(OR) + CI2 k = 


= 0.002662 


ooy ) 


OLRl(OR) + Cro2 OLRl(CR) k 


= 0.314159 




OLRl(CR) ^ OLRl(OR) + Cro2 k 


= 0.068319 


oyi ) 


OLRl(OC) + CI2 ^ OLRl(RC) k = 


= 0.314159 


392) 


OLRl(RC) ^ OLRl(OC) + CI2 k = 


= 0.346100 


393) 


OLRl(OC) + Cro2 ^ OLRl(CC) k 


= 0.314159 


394) 


OLRl(CC) ^ OLRl(OC) + Croa k 


= 0.068319 


395) 


OLRl(ORp) + CI2 ^ OLRl(RRp) k = 0.314159 
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396) OLRl(RRp) ^ OLRl(0Rp)+ CI2 k = 0.346100 

397) OLRl(ORp) + Cro2 ^ OLRl(CRp) k = 0.314159 

398) OLRl(CRp) ^ OLRl(ORp) + Cro2 k = 0.068319 

399) OLRl(ORp) OLRl(OO) + Rp + MCI k = 0.011000 

400) OLRl(RRp) ^ OLRl(RO) + Rp + MCI k = 0.011000 

401) OLRl(CRp) ^ OLRl(CO) + Rp + MCI k = 0.011000 

402) OL(RRO) + O(ORR) ^ OLR2(00) k = 62.095095 

403) OLR2(00) ^ OL(RRO) + O(ORR) k = 0.008046 

404) OL(RRO) + O(RRR) ^ OLR2(0R) k = 62.095095 

405) OLR2(0R) ^ OL(RRO) + O(RRR) k = 0.008046 

406) OL(RRO) + O(CRR) ^ OLR2(0C) k = 62.095095 

407) OLR2(0C) ^ OL(RRO) + O(CRR) k = 0.008046 

408) OL(RRO) + O(RpRR) ^ OLR2(0Rp) k = 62.095095 

409) OLR2(0Rp) ^ OL(RRO) + O(RpRR) k = 0.008046 

410) OL(RRR) + O(ORR) ^ OLR2(R0) k = 62.095095 

411) OLR2(R0) ^ OL(RRR) + O(ORR) k = 0.008046 

412) OL(RRR) + O(RRR) ^ 0LR2(RR) k = 62.095095 

413) 0LR2(RR) ^ OL(RRR) + O(RRR) k = 0.008046 

414) OL(RRR) + O(CRR) ^ 0LR2(RC) k = 62.095095 

415) 0LR2(RC) ^ OL(RRR) + O(CRR) k = 0.008046 

416) OL(RRR) + O(RpRR) ^ 0LR2(RRp) k = 62.095095 

417) 0LR2(RRp) ^ OL(RRR) + O(RpRR) k = 0.008046 

418) OL(RRC) + O(ORR) ^ OLR2(C0) k = 62.095095 

419) OLR2(C0) ^ OL(RRC) + O(ORR) k = 0.008046 

420) OL(RRC) + O(RRR) ^ 0LR2(CR) k = 62.095095 

421) 0LR2(CR) ^ OL(RRC) + O(RRR) k = 0.008046 

422) OL(RRC) + O(CRR) ^ 0LR2(CC) k = 62.095095 

423) 0LR2(CC) ^ OL(RRC) + O(CRR) k = 0.008046 

424) OL(RRC) + O(RpRR) ^ 0LR2(CRp) k = 62.095095 

425) 0LR2(CRp) ^ OL(RRC) + O(RpRR) k = 0.008046 

426) OL(RRRp) + O(ORR) ^ OLR2(RpO) k = 62.095095 

427) OLR2(RpO) ^ OL(RRRp) + O(ORR) k = 0.008046 

428) OL(RRRp) + O(RRR) ^ 0LR2(RpR) k = 62.095095 
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A o^^ 
4zy J 


Uijriz(^±tprLj — > Uij(^ri,rtxlpj + U^^xtrtxlj K — U.UUoU4D 


A Q^^ 
4oU J 


vJijt^xtKrlpj + U^^^xlKj — > U-Lxtz^KpUj K — Dz.UyOUyO 


/I Q 1 ^ 
4oi J 


0LR2(RpC) ^ OL(RRRp) + O(CRR) k = 0.008046 


4oz J 


OL(RRRp) + O(RpRR) ^ 0LR2(RpRp) k = 62.095095 


AQQ\ 

4oc) J 


0LR2(RpRp) ^ OL(RRRp) + O(RpRR) k = 0.008046 


4o4 J 


OLR2(00) + CI2 OLR2(R0) k = 


0.314159 


4ob j 


OLR2(R0) OLR2(00) + CI2 k = 


0.035701 


4oD j 


OLR2(00) + CI2 OLR2(0R) k = 


0.314159 


4o r j 


OLR2(0R) OLR2(00) + CI2 k = 


38.256936 


4oo j 


OLR2(00) + Cro2 OLR2(C0) k = 


= 0.314159 


4oy J 


OLR2(C0) OLR2(00) + Cro2 k = 


= 0.662315 


44U J 


OLR2(00) + Cro2 OLR2(0C) k = 


= 0.314159 


44i J 


OLR2(0C) OLR2(00) + Cro2 k = 


= 0.068319 


44z J 


OLR2(00) + Rp ^ OLR2(0Rp) k = 


= 0.314159 


44o J 


OLR2(0Rp) OLR2(00) + Rp k = 


= 1.490712 


/I /I /I ^ 
444 J 


OLR2(00) + Rp ^ OLR2(RpO) k = 


= 0.314159 


44b J 


OLR2(RpO) OLR2(00) + Rp k = 


= 2.062175 


44D J 


OLR2(R0) + CI2 0LR2(RR) k = 


= 0.314159 


44 ( J 


0LR2(RR) OLR2(R0) + CI2 k = 


= 38.256936 


/I /I s^ 

448 J 


OLR2(R0) + Cro2 0LR2(RC) k 


= 0.314159 


44y J 


0LR2(RC) OLR2(R0) + Cro2 k 


= 0.068319 


4oU ) 


OLR2(R0) + Rp ^ 0LR2(RRp) k 


= 0.314159 


AK^ \ 
4oi J 


0LR2(RRp) OLR2(R0) + Rp k 


= 1.490712 


4oz J 


OLR2(C0) + CI2 0LR2(CR) k = 


= 0.314159 


4oo j 


UHxZ[(^ix) — > UHxZ[y^{)) + Ui2 K = 


= OO.ZOoyoD 


404 J 


OLR2(C0) + Cro2 ^ 0LR2(CC) k 


= 0.314159 


4ot) j 


0LR2(CC) ^ OLR2(C0) + Cros k 


= 0.068319 


4dd j 


OLR2(C0) + Rp ^ 0LR2(CRp) k 


= 0.314159 


40 


0LR2(CRp) ^ OLR2(C0) + Rp k 


= 1.490712 


458) 


OLR2(0R) + CI2 ^ 0LR2(RR) k = 


= 0.314159 


459) 


0LR2(RR) ^ OLR2(0R) + CI2 k = 


= 0.035701 


460) 


OLR2(0R) + Cro2 ^ 0LR2(CR) k 


= 0.314159 


461) 


0LR2(CR) ^ OLR2(0R) + Cro2 k 


= 0.662315 



462) OLR2(0R) + Rp ^ 0LR2(RpR) k = 0.314159 

463) 0LR2(RpR) ^ OLR2(0R) + Rp k = 2.062175 

464) OLR2(0C) + CI2 ^ 0LR2(RC) k = 0.314159 

465) 0LR2(RC) ^ OLR2(0C) + CI2 k = 0.035701 

466) OLR2(0C) + Croa ^ 0LR2(CC) k = 0.314159 

467) 0LR2(CC) ^ OLR2(0C) + Cros k = 0.662315 

468) OLR2(0C) + Rp ^ 0LR2(RpC) k = 0.314159 

469) 0LR2(RpC) ^ OLR2(0C) + Rp k = 2.062175 

470) OLR2(0Rp) + CI2 ^ 0LR2(RRp) k = 0.314159 

471) 0LR2(RRp) ^ OLR2(0Rp) + CI2 k = 0.035701 

472) OLR2(0Rp) + Cro2 ^ 0LR2(CRp) k = 0.314159 

473) 0LR2(CRp) OLR2(0Rp) + Cro2 k = 0.662315 

474) OLR2(0Rp) + Rp ^ 0LR2(RpRp) k = 0.314159 

475) 0LR2(RpRp) ^ OLR2(0Rp) + Rp k = 2.062175 

476) OLR2(0Rp) ^ OLR2(00) + Rp + MCI k = 0.011000 

477) 0LR2(RRp) ^ OLR2(R0) + Rp + MCI k = 0.011000 

478) 0LR2(CRp) ^ OLR2(C0) + Rp + MCI k = 0.011000 

479) 0LR2(RpRp) ^ OLR2(RpO) + Rp + MCI k = 0.011000 

480) OL(ORR) + 0(RRO ^ OLR3(00) k = 62.095095 

481) OLR3(00) ^ OL(ORR) + O(RRO) k = 0.008046 

482) OL(0RR + 0(RRR ^ OLR3(0R) k = 62.095095 

483) OLR3(0R) ^ OL(ORR) + O(RRR) k = 0.008046 

484) OL(ORR) + 0(RRC ^ OLR3(0C) k = 62.095095 

485) OLR3(0C) ^ OL(ORR) + O(RRC) k = 0.008046 

486) OL(RRR) + 0(RRO OLR3(R0) k = 62.095095 

487) OLR3(R0) ^ OL(RRR) + O(RRO) k = 0.008046 

488) OL(RRR) + 0(RRR ^ 0LR3(RR) k = 62.095095 

489) 0LR3(RR) ^ OL(RRR) + O(RRR) k = 0.008046 

490) OL(RRR) + 0(RRC ^ 0LR3(RC) k = 62.095095 

491) 0LR3(RC) ^ OL(RRR) + O(RRC) k = 0.008046 

492) OL(CRR) + 0(RRO ^ OLR3(C0) k = 62.095095 

493) OLR3(C0) ^ OL(CRR) + O(RRO) k = 0.008046 

494) OL(CRR) + 0(RRR ^ 0LR3(CR) k = 62.095095 



495) 0LR3(CR) ^ OL(CRR) + O(RRR) k = 0.008046 

496) OL(CRR) + 0(RRC ^ 0LR3(CC) k = 62.095095 

497) 0LR3(CC) ^ OL(CRR) + O(RRC) k = 0.008046 

498) OLR3(00) + CI2 ^ OLR3(R0) k = 0.314159 

499) OLR3(R0) ^ OLR3(00) + CI2 k = 0.346100 

500) OLR3(00) + CI2 ^ OLR3(0R) k = 0.314159 

501) OLR3(0R) ^ OLR3(00) + CI2 k = 0.294263 

502) OLR3(00) + Cro2 ^ OLR3(C0) k = 0.314159 

503) OLR3(C0) ^ OLR3(00) + Cro2 k = 0.068319 

504) OLR3(00) + Cro2 ^ OLR3(0C) k = 0.314159 

505) OLR3(0C) ^ OLR3(00) + Cro2 k = 0.662315 

506) OLR3(R0) + CI2 ^ 0LR3(RR) k = 0.314159 

507) 0LR3(RR) ^ OLR3(R0) + CI2 k = 0.294263 

508) OLR3(R0) + Cros ^ 0LR3(RC) k = 0.314159 

509) 0LR3(RC) ^ OLR3(R0) + Cros k = 0.662315 

510) OLR3(C0) + CI2 ^ 0LR3(CR) k = 0.314159 

511) 0LR3(CR) ^ OLR3(C0) + CI2 k = 0.294263 

512) OLR3(C0) + Cro2 ^ 0LR3(CC) k = 0.314159 

513) 0LR3(CC) ^ OLR3(C0) + Cro2 k = 0.662315 

514) OLR3(0R) + CI2 ^ 0LR3(RR) k = 0.314159 

515) 0LR3(RR) ^ OLR3(0R) + CI2 k = 0.346100 

516) OLR3(0R) + Cro2 ^ 0LR3(CR) k = 0.314159 

517) 0LR3(CR) ^ OLR3(0R) + Cro2 k = 0.068319 

518) OLR3(0C) + CI2 ^ 0LR3(RC) k = 0.314159 

519) 0LR3(RC) ^ OLR3(0C) + CI2 k = 0.346100 

520) OLR3(0C) + Cro2 ^ 0LR3(CC) k = 0.314159 

521) 0LR3(CC) ^ OLR3(0C) + Cro2 k = 0.068319 

522) OL(RRO) + 0(RRO ^ OLR4(00) k = 62.095095 

523) OLR4(00) ^ OL(RRO) + O(RRO) k = 0.008046 

524) OL(RRO) + 0(RRR ^ OLR4(0R) k = 62.095095 

525) OLR4(0R) ^ OL(RRO) + O(RRR) k = 0.008046 

526) OL(RRO) + 0(RRC ^ OLR4(0C) k = 62.095095 

527) OLR4(0C) ^ OL(RRO) + O(RRC) k = 0.008046 
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528) OL(RRR) + 0(RRO ^ OLR4(R0) k = 62.095095 

529) OLR4(R0) ^ OL(RRR) + O(RRO) k = 0.008046 

530) OL(RRR) + 0(RRR ^ 0LR4(RR) k = 62.095095 

531) 0LR4(RR) ^ OL(RRR) + O(RRR) k = 0.008046 

532) OL(RRR) + 0(RRC ^ 0LR4(RC) k = 62.095095 

533) 0LR4(RC) ^ OL(RRR) + O(RRC) k = 0.008046 

534) OL(RRC) + 0(RRO ^ OLR4(C0) k = 62.095095 

535) OLR4(C0) ^ OL(RRC) + O(RRO) k = 0.008046 

536) OL(RRC) + 0(RRR ^ 0LR4(CR) k = 62.095095 

537) 0LR4(CR) ^ OL(RRC) + O(RRR) k = 0.008046 

538) OL(RRC) + 0(RRC ^ 0LR4(CC) k = 62.095095 

539) 0LR4(CC) ^ OL(RRC) + O(RRC) k = 0.008046 

540) OL(RRRp) + 0(RRO ^ OLR4(RpO) k = 62.095095 

541) OLR4(RpO) ^ OL(RRRp) + O(RRO) k = 0.008046 

542) OL(RRRp) + 0(RRR ^ 0LR4(RpR) k = 62.095095 

543) 0LR4(RpR) ^ OL(RRRp) + O(RRR) k = 0.008046 

544) OL(RRRp) + 0(RRC ^ 0LR4(RpC) k = 62.095095 

545) 0LR4(RpC ^ OL(RRRp) + O(RRC) k = 0.008046 

546) OLR4(00) + CI2 ^ OLR4(R0) k = 0.314159 

547) OLR4(R0) ^ OLR4(00) + CI2 k = 0.035701 

548) OLR4(00) + CI2 ^ OLR4(0R) k = 0.314159 

549) OLR4(0R) ^ OLR4(00) + CI2 k = 0.294263 

550) OLR4(00) + Cro2 ^ OLR4(C0) k = 0.314159 

551) OLR4(C0) ^ OLR4(00) + Cro2 k = 0.662315 

552) OLR4(00) + Cro2 ^ OLR4(0C) k = 0.314159 

553) OLR4(0C) ^ OLR4(00) + Cro2 k = 0.662315 

554) OLR4(00) + Rp ^ OLR4(RpO) k = 0.314159 

555) OLR4(RpO) ^ OLR4(00) + Rp k = 2.062175 

556) OLR4(R0) + CI2 ^ 0LR4(RR) k = 0.314159 

557) 0LR4(RR) ^ OLR4(R0) + CI2 k = 0.002263 

558) OLR4(R0) + Cro2 ^ 0LR4(RC) k = 0.314159 

559) 0LR4(RC) ^ OLR4(R0) + Cro2 k = 0.662315 

560) OLR4(C0) + CI2 ^ 0LR4(CR) k = 0.314159 
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oui ) 


OT T5 A fr^T) \ 

Uljri4(^Uri,j — 


i r\TT>A(nn\ _l ht - 


- n oo/iofiQ 

- U.zy4z0o 


00/ j 




Ur02 — > ULiix^[y-iKj ) K 


— U.oi4ioy 


ODO ) 


OT T5 /I tr~'(~'\ 

UijK4(UL'j — 


> UijK4^UUj + Ur02 K 


= U.oozolo 


004 ) 


<JijK4(^UKj + 


r^j , OT T3 /I /T3 T3 ^ 1^ — 

Ui2 — > KjHx^[iXSX) K - 


- U.oi4ioy 


000 ) 




> W-LK4(^UKj + Ui2 K - 


- U.UUUZ/O 


OOD ) 


UljK4(UKj + 


Ur02 — > KjHx4i[LyIx) K 


n Q 1 /1 1 F;n 
= U.oi4ioy 


00 ( j 




> Uijxt4(^Uxtj + (-jT02 K 


— U.oozoio 


OOo ) 


Uijri4(^UKj + 


rtp — > Uljxl4(^xlprtj K 


— U.oi4ioy 


ooy j 




— UijK4(UKj + xlp K : 


= z.Uozi 


/Uj 


Uijri4(^UU ) + 




- n QI /1 1 KCi 

- u.oi4ioy 


/ i j 


Uijri4(^rlU j — 


> KjLiX^[(j\^ ) + Ui2 K - 


- n nQf^vni 

- U.UoO /Ul 


572) 


OLR4(0C) + 


Cro2 0LR4(CC) k 


= 0.314159 


573) 


OLR4(0C) ^ 


0LR4(CC) + Cro2 k 


= 0.662315 


574) 


OLR4(0C) + 


Rp ^ 0LR4(RpC) k : 


= 0.314159 


575) 


0LR4(RpC) - 


OLR4(0C) + Rp k : 


= 2.062175 



Including non-specific binding extends the reaction to scheme with these four reactions (rate 
constants calculated for AGnsb = — 3.5 kcal/mol): 

576) D + CI2 ^ 1 DCI k = 0.314159 

577) DCI ^ 1 D + 1 CI2 k = 646634.937870 

578) D + Cro2 1 DCRO k = 0.314159 

579) DCRO -> 1 D + 1 Cro2 k = 646634.937870 
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